LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
associationUtils.py
Go to the documentation of this file.
1# This file is part of pipe_tasks.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21
22"""Utilities for interfacing with hpgeom. Originally implemented in
23http://github.com/LSSTDESC/dia_pipe and then translated to hpgeom.
24"""
25
26__all__ = ["toIndex", "toRaDec", "eq2xyz", "eq2xyzVec", "convert_spherical",
27 "convert_spherical_array", "query_disc"]
28
29import hpgeom as hpg
30import numpy as np
31
32
33def toIndex(nside, ra, dec):
34 """Return healpix index given ra, dec in degrees
35
36 Parameters
37 ----------
38 nside : `int`
39 Power of 2 nside healpix resolution.
40 ra : `float`
41 RA in degrees.
42 dec : `float`
43 Declination in degrees
44
45 Returns
46 -------
47 index : `int`
48 Unique healpix pixel ID containing point RA, DEC at resolution nside.
49 """
50 return hpg.angle_to_pixel(nside, ra, dec, nest=False)
51
52
53def toRaDec(nside, index):
54 """Convert from healpix index to ra,dec in degrees
55
56 Parameters
57 ----------
58 nside : `int`
59 Resolution of healpixel "grid".
60 index : `int`
61 Index of the healpix pixel we want to find the location of.
62
63 Returns
64 -------
65 pos : `numpy.ndarray`, (2,)
66 RA and DEC of healpix pixel location in degrees.
67 """
68 ra, dec = hpg.pixel_to_angle(nside, index, nest=False)
69 return np.dstack((ra, dec))[0]
70
71
72def eq2xyz(ra, dec):
73 """Convert from equatorial ra,dec in degrees to x,y,z on unit sphere.
74
75 Parameters
76 ----------
77 ra : `float`
78 RA in degrees.
79 dec : `float`
80 Declination in degrees
81
82 Returns
83 -------
84 xyz : `numpy.ndarray`, (3,)
85 Float xyz positions on the unit sphere.
86 """
87 phi = np.deg2rad(ra)
88 theta = np.pi/2 - np.deg2rad(dec)
89 sintheta = np.sin(theta)
90 x = sintheta*np.cos(phi)
91 y = sintheta*np.sin(phi)
92 z = np.cos(theta)
93 return np.array([x, y, z])
94
95
96def eq2xyzVec(ra, dec):
97 """Convert equatorial ra,dec in degrees to x,y,z on the unit sphere
98 parameters
99
100 Vectorized version of ``eq2xyz``
101
102 Parameters
103 ----------
104 ra : array_like, (N,)
105 Array of RA in degrees.
106 dec : array_like, (N,)
107 Declination in degrees
108
109 Returns
110 -------
111 vec : `numpy.ndarray`, (N,3)
112 Array of unitsphere 3-vectors.
113 """
114 ra = np.array(ra, dtype='f8', ndmin=1, copy=False)
115 dec = np.array(dec, dtype='f8', ndmin=1, copy=False)
116 if ra.size != dec.size:
117 raise ValueError("ra,dec not same size: %s,%s" % (ra.size, dec.size))
118
119 vec = eq2xyz(ra, dec)
120
121 return vec
122
123
124def convert_spherical(ra, dec):
125 """Convert from ra,dec to spherical coordinates.
126
127 Used in query_disc.
128
129 Parameters
130 ----------
131 ra : `float`
132 RA in radians.
133 dec : `float`
134 Declination in radians
135 """
136 return np.dstack([np.cos(dec*np.pi/180)*np.cos(ra*np.pi/180),
137 np.cos(dec*np.pi/180)*np.sin(ra*np.pi/180),
138 np.sin(dec*np.pi/180)])[0]
139
140
142 """Convert from and a array ra,dec to spherical coordinates.
143
144 Used in query_disc
145
146 Parameters
147 ----------
148 array : `numpy.ndarray`, (N, 2)
149 (N, 2) Array of RA, DEC values.
150
151 Returns
152 -------
153 vecs : `numpy.ndarray`, (N, 3)
154 Vectors on the unit sphere
155 """
156 ra = array[:, 0]
157 dec = array[:, 1]
158 return convert_spherical(ra, dec)
159
160
161def query_disc(nside, ra, dec, max_rad, min_rad=0):
162 """Get the list of healpix indices within max_rad, min_rad given in radians
163 around ra,dec given in degrees
164
165 Parameters
166 ----------
167 nside : `int`
168 Resolution of the healpixels to search/return.
169 ra : `float`
170 RA in degrees.
171 dec : `float`
172 Declination in degrees
173 max_rad : `float`
174 Max distance in radians to search nearby healpixels.
175 min_rad : `float`, optional
176 Minimum distance in radians to search healpixels. Default = 0.
177 """
178 ra = np.atleast_1d(ra)
179 dec = np.atleast_1d(dec)
180
181 max_rad_deg = np.rad2deg(max_rad)
182
183 pixels = np.unique(
184 [hpg.query_circle(nside, a, b, max_rad_deg, nest=False)
185 for (a, b) in zip(ra, dec)])
186
187 if min_rad > 0 and len(pixels) > 0:
188 vec0 = convert_spherical(ra, dec)
189 min_rad2 = min_rad**2
190 vecs = convert_spherical_array(toRaDec(nside, pixels))
191 dsq = np.sum((vecs - vec0)**2, axis=1)
192 match = dsq > min_rad2
193 pixels = pixels[match]
194
195 return pixels
query_disc(nside, ra, dec, max_rad, min_rad=0)