Coverage for python/lsst/pipe/tasks/associationUtils.py: 9%

1# This file is part of pipe_tasks.

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.

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/>.

22"""Utilities for interfacing with hpgeom. Originally implemented in

23http://github.com/LSSTDESC/dia_pipe and then translated to hpgeom.

24"""

26__all__ = ["toIndex", "toRaDec", "eq2xyz", "eq2xyzVec", "convert_spherical",

27 "convert_spherical_array", "query_disc", "obj_id_to_ss_object_id", "ss_object_id_to_obj_id"]

29from astropy.time import Time

30import hpgeom as hpg

31import numpy as np

32import numpy.typing as npt

35def toIndex(nside, ra, dec):

36 """Return healpix index given ra, dec in degrees

38 Parameters

39 ----------

40 nside : `int`

41 Power of 2 nside healpix resolution.

42 ra : `float`

43 RA in degrees.

44 dec : `float`

45 Declination in degrees

47 Returns

48 -------

49 index : `int`

50 Unique healpix pixel ID containing point RA, DEC at resolution nside.

51 """

52 return hpg.angle_to_pixel(nside, ra, dec, nest=False)

55def toRaDec(nside, index):

56 """Convert from healpix index to ra,dec in degrees

58 Parameters

59 ----------

60 nside : `int`

61 Resolution of healpixel "grid".

62 index : `int`

63 Index of the healpix pixel we want to find the location of.

65 Returns

66 -------

67 pos : `numpy.ndarray`, (2,)

68 RA and DEC of healpix pixel location in degrees.

69 """

70 ra, dec = hpg.pixel_to_angle(nside, index, nest=False)

71 return np.dstack((ra, dec))[0]

74def eq2xyz(ra, dec):

75 """Convert from equatorial ra,dec in degrees to x,y,z on unit sphere.

77 Parameters

78 ----------

79 ra : `float`

80 RA in degrees.

81 dec : `float`

82 Declination in degrees

84 Returns

85 -------

86 xyz : `numpy.ndarray`, (3,)

87 Float xyz positions on the unit sphere.

88 """

89 phi = np.deg2rad(ra)

90 theta = np.pi/2 - np.deg2rad(dec)

91 sintheta = np.sin(theta)

92 x = sintheta*np.cos(phi)

93 y = sintheta*np.sin(phi)

94 z = np.cos(theta)

95 return np.array([x, y, z])

98def eq2xyzVec(ra, dec):

99 """Convert equatorial ra,dec in degrees to x,y,z on the unit sphere

100 parameters

101

102 Vectorized version of ``eq2xyz``

103

104 Parameters

105 ----------

106 ra : array_like, (N,)

107 Array of RA in degrees.

108 dec : array_like, (N,)

109 Declination in degrees

110

111 Returns

112 -------

113 vec : `numpy.ndarray`, (N,3)

114 Array of unitsphere 3-vectors.

115 """

116 ra = np.array(ra, dtype='f8', ndmin=1, copy=False)

117 dec = np.array(dec, dtype='f8', ndmin=1, copy=False)

118 if ra.size != dec.size:

119 raise ValueError("ra,dec not same size: %s,%s" % (ra.size, dec.size))

120

121 vec = eq2xyz(ra, dec)

122

123 return vec

124

125

126def convert_spherical(ra, dec):

127 """Convert from ra,dec to spherical coordinates.

128

129 Used in query_disc.

130

131 Parameters

132 ----------

133 ra : `float`

134 RA in radians.

135 dec : `float`

136 Declination in radians

137 """

138 return np.dstack([np.cos(dec*np.pi/180)*np.cos(ra*np.pi/180),

139 np.cos(dec*np.pi/180)*np.sin(ra*np.pi/180),

140 np.sin(dec*np.pi/180)])[0]

141

142

143def convert_spherical_array(array):

144 """Convert from and a array ra,dec to spherical coordinates.

145

146 Used in query_disc

147

148 Parameters

149 ----------

150 array : `numpy.ndarray`, (N, 2)

151 (N, 2) Array of RA, DEC values.

152

153 Returns

154 -------

155 vecs : `numpy.ndarray`, (N, 3)

156 Vectors on the unit sphere

157 """

158 ra = array[:, 0]

159 dec = array[:, 1]

160 return convert_spherical(ra, dec)

161

162

163def query_disc(nside, ra, dec, max_rad, min_rad=0):

164 """Get the list of healpix indices within max_rad, min_rad given in radians

165 around ra,dec given in degrees

166

167 Parameters

168 ----------

169 nside : `int`

170 Resolution of the healpixels to search/return.

171 ra : `float`

172 RA in degrees.

173 dec : `float`

174 Declination in degrees

175 max_rad : `float`

176 Max distance in radians to search nearby healpixels.

177 min_rad : `float`, optional

178 Minimum distance in radians to search healpixels. Default = 0.

179 """

180 ra = np.atleast_1d(ra)

181 dec = np.atleast_1d(dec)

182

183 max_rad_deg = np.rad2deg(max_rad)

184

185 pixels = np.unique(

186 [hpg.query_circle(nside, a, b, max_rad_deg, nest=False)

187 for (a, b) in zip(ra, dec)])

188

189 if min_rad > 0 and len(pixels) > 0:

190 vec0 = convert_spherical(ra, dec)

191 min_rad2 = min_rad**2

192 vecs = convert_spherical_array(toRaDec(nside, pixels))

193 dsq = np.sum((vecs - vec0)**2, axis=1)

194 match = dsq > min_rad2

195 pixels = pixels[match]

196

197 return pixels

198

199

200def obj_id_to_ss_object_id(objID):

201 """Convert from Minor Planet Center packed provisional object ID to

202 Rubin ssObjectID.

203

204 Parameters

205 ----------

206 objID : `str`

207 Minor Planet Center packed provisional designation for a small solar

208 system object. Must be fewer than eight characters.

209

210 Returns

211 -------

212 ssObjectID : `int`

213 Rubin ssObjectID

214

215 Raises

216 ------

217 ValueError

218 Raised if either objID is shorter than 7 or longer than 8 characters or contains

219 illegal objID characters

220 """

221 if len(objID) > 8:

222 raise ValueError(f'objID longer than 8 characters: "{objID}"')

223 if len(objID) < 7:

224 raise ValueError(f'objID shorter than 7 characters: "{objID}"')

225 if any([ord(c) > 255 for c in objID]):

226 raise ValueError(f'{[c for c in objID if ord(c) > 255]} not legal objID characters (ascii [1, 255])')

227

228 ssObjectID = ord(objID[0])

229 for character in objID[1:]:

230 ssObjectID <<= 8

231 ssObjectID += ord(character)

232 return ssObjectID

233

234

235def ss_object_id_to_obj_id(ssObjectID):

236 """Convert from Rubin ssObjectID to Minor Planet Center packed provisional

237 object ID.

238

239 Parameters

240 ----------

241 ssObjectID : `int`

242 Rubin ssObjectID

243

244 Returns

245 -------

246 objID : `str`

247 Minor Planet Center packed provisional designation.

248

249 Raises

250 ------

251 """

252 if ssObjectID < 0 or ssObjectID >= (1 << 64):

253 raise ValueError(f'ssObjectID ({ssObjectID}) outside [0, 2^64 - 1].')

254

255 objID = ''.join([chr((ssObjectID >> (8 * i)) % 256) for i in reversed(range(0, 8))])

256 objID = objID.replace('\x00', '')

257 return objID

258

259# All the below designation-related code are copied from B612's adam_core

260# adam_core should eventually be added as an external dependency, and this

261# should be replaced with imports.

262

263

264def pack_mpc_designation(designation: str) -> str:

265 """

266 B612 code to pack a unpacked MPC designation. For example,

267 provisional designation 1998 SS162 will be packed to J98SG2S.

268 Permanent designation 323 will be packed to 00323.

269

270 TODO: add support for comet and natural satellite designations

271

272 Parameters

273 ----------

274 designation : str

275 MPC unpacked designation.

276

277 Returns

278 -------

279 designation_pf : str

280 MPC packed form designation.

281

282 Raises

283 ------

284 ValueError : If designation cannot be packed.

285 """

286 # Lets see if its a numbered object

287 try:

288 return pack_numbered_designation(designation)

289 except ValueError:

290 pass

291

292 # If its not numbered, maybe its a provisional designation

293 try:

294 return pack_provisional_designation(designation)

295 except ValueError:

296 pass

297

298 # If its a survey designation, deal with it

299 try:

300 return pack_survey_designation(designation)

301 except ValueError:

302 pass

303

304 err = (

305 "Unpacked designation '{}' could not be packed.\n"

306 "It could not be recognized as any of the following:\n"

307 " - a numbered object (e.g. '3202', '203289', '3140113')\n"

308 " - a provisional designation (e.g. '1998 SV127', '2008 AA360')\n"

309 " - a survey designation (e.g. '2040 P-L', '3138 T-1')"

310 )

311 raise ValueError(err.format(designation))

312

313

314BASE62 = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"

315BASE62_MAP = {BASE62[i]: i for i in range(len(BASE62))}

316

317

318def _unpack_mpc_date(epoch_pf: str) -> Time:

319 # Taken from Lynne Jones' SSO TOOLS.

320 # See https://minorplanetcenter.net/iau/info/PackedDates.html

321 # for MPC documentation on packed dates.

322 # Examples:

323 # 1998 Jan. 18.73 = J981I73

324 # 2001 Oct. 22.138303 = K01AM138303

325 epoch_pf = str(epoch_pf)

326 year = int(epoch_pf[0], base=32) * 100 + int(epoch_pf[1:3])

327 month = int(epoch_pf[3], base=32)

328 day = int(epoch_pf[4], base=32)

329 isot_string = "{:d}-{:02d}-{:02d}".format(year, month, day)

330

331 if len(epoch_pf) > 5:

332 fractional_day = float("." + epoch_pf[5:])

333 hours = int((24 * fractional_day))

334 minutes = int(60 * ((24 * fractional_day) - hours))

335 seconds = 3600 * (24 * fractional_day - hours - minutes / 60)

336 isot_string += "T{:02d}:{:02d}:{:09.6f}".format(hours, minutes, seconds)

337

338 return Time(isot_string, format="isot", scale="tt")

339

340

341def convert_mpc_packed_dates(pf_tt: npt.ArrayLike) -> Time:

342 """

343 Convert MPC packed form dates (in the TT time scale) to

344 MJDs in TT. See: https://minorplanetcenter.net/iau/info/PackedDates.html

345 for details on the packed date format.

346

347 Parameters

348 ----------

349 pf_tt : `~numpy.ndarray` (N)

350 MPC-style packed form epochs in the TT time scale.

351

352 Returns

353 -------

354 mjd_tt : `~astropy.time.core.Time` (N)

355 Epochs in TT MJDs.

356 """

357 isot_tt = []

358 for epoch in pf_tt:

359 isot_tt.append(_unpack_mpc_date(epoch))

360

361 return Time(isot_tt)

362

363

364def pack_numbered_designation(designation: str) -> str:

365 """

366 Pack a numbered MPC designation.

367

368 Examples of numbered designations:

369 Numbered Packed

370 3202 03202

371 50000 50000

372 100345 A0345

373 360017 a0017

374 203289 K3289

375 620000 ~0000

376 620061 ~000z

377 3140113 ~AZaz

378 15396335 ~zzzz

379

380 Parameters

381 ----------

382 designation : str

383 MPC numbered designation.

384

385 Returns

386 -------

387 designation_pf : str

388 MPC packed numbered designation.

389

390 Raises

391 ------

392 ValueError : If the numbered designation cannot be packed.

393 If the numbered designation is larger than 15396335.

394 """

395 number = int(designation)

396 if number > 15396335:

397 raise ValueError(

398 "Numbered designation is too large. Maximum supported is 15396335."

399 )

400

401 if number <= 99999:

402 return "{:05}".format(number)

403 elif (number >= 100000) and (number <= 619999):

404 bigpart, remainder = divmod(number, 10000)

405 return f"{BASE62[bigpart]}{remainder:04}"

406 else:

407 x = number - 620000

408 number_pf = []

409 while x:

410 number_pf.append(BASE62[int(x % 62)])

411 x //= 62

412

413 number_pf.reverse()

414 return "~{}".format("".join(number_pf).zfill(4))

415

416

417def pack_provisional_designation(designation: str) -> str:

418 """

419 Pack a provisional MPC designation.

420

421 Examples of provisional designations:

422 Provisional Packed

423 1995 XA J95X00A

424 1995 XL1 J95X01L

425 1995 FB13 J95F13B

426 1998 SQ108 J98SA8Q

427 1998 SV127 J98SC7V

428 1998 SS162 J98SG2S

429 2099 AZ193 K99AJ3Z

430 2008 AA360 K08Aa0A

431 2007 TA418 K07Tf8A

432

433 Parameters

434 ----------

435 designation : str

436 MPC provisional designation.

437

438 Returns

439 -------

440 designation_pf : str

441 MPC packed provisional designation.

442

443 Raises

444 ------

445 ValueError : If the provisional designation cannot be packed.

446 The provisional designations is not at least 6 characters long.

447 The first 4 characters of the provisional designation are not a year.

448 The 5th character of the provisional designation is not a space.

449 The provisional designation contains a hyphen.

450 The half-month letter is I or Z.

451 """

452 if len(designation) < 6:

453 raise ValueError(

454 "Provisional designations should be at least 6 characters long."

455 )

456 if not designation[:3].isdecimal():

457 raise ValueError(

458 "Expected the first 4 characters of the provisional designation to be a year."

459 )

460 if designation[4] != " ":

461 raise ValueError(

462 "Expected the 5th character of the provisional designation to be a space."

463 )

464 if "-" in designation:

465 raise ValueError("Provisional designations cannot contain a hyphen.")

466

467 year = BASE62[int(designation[0:2])] + designation[2:4]

468 letter1 = designation[5]

469 letter2 = designation[6]

470 cycle = designation[7:]

471

472 if letter1 in {"I", "Z"}:

473 raise ValueError("Half-month letters cannot be I or Z.")

474 if letter1.isdecimal() or letter2.isdecimal():

475 raise ValueError("Invalid provisional designation.")

476

477 cycle_pf = "00"

478 if len(cycle) > 0:

479 cycle_int = int(cycle)

480 if cycle_int <= 99:

481 cycle_pf = str(cycle_int).zfill(2)

482 else:

483 cycle_pf = BASE62[cycle_int // 10] + str(cycle_int % 10)

484

485 designation_pf = "{}{}{}{}".format(year, letter1, cycle_pf, letter2)

486 return designation_pf

487

488

489def pack_survey_designation(designation: str) -> str:

490 """

491 Pack a survey MPC designation.

492

493 Examples of survey designations:

494 Survey Packed

495 2040 P-L PLS2040

496 3138 T-1 T1S3138

497 1010 T-2 T2S1010

498 4101 T-3 T3S4101

499

500 Parameters

501 ----------

502 designation : str

503 MPC survey designation.

504

505 Returns

506 -------

507 designation_pf : str

508 MPC packed survey designation.

509

510 Raises

511 ------

512 ValueError : If the survey designation cannot be packed.

513 The survey designation does not start with P-L, T-1, T-2, or T-3.

514 """

515 number = designation[0:4]

516 survey = designation[5:]

517

518 if survey == "P-L":

519 survey_pf = "PLS"

520

521 elif survey[0:2] == "T-" and survey[2] in {"1", "2", "3"}:

522 survey_pf = "T{}S".format(survey[2])

523

524 else:

525 raise ValueError("Survey designations must start with P-L, T-1, T-2, T-3.")

526

527 designation_pf = "{}{}".format(survey_pf, number.zfill(4))

528 return designation_pf

529

530

531def unpack_numbered_designation(designation_pf: str) -> str:

532 """

533 Unpack a numbered MPC designation.

534

535 Examples of numbered designations:

536 Numbered Unpacked

537 03202 3202

538 50000 50000

539 A0345 100345

540 a0017 360017

541 K3289 203289

542 ~0000 620000

543 ~000z 620061

544 ~AZaz 3140113

545 ~zzzz 15396335

546

547 Parameters

548 ----------

549 designation_pf : str

550 MPC packed numbered designation.

551

552 Returns

553 -------

554 designation : str

555 MPC unpacked numbered designation.

556

557 Raises

558 ------

559 ValueError : If the numbered designation cannot be unpacked.

560 The packed numbered designation is not at least 4 characters long.

561 """

562 number = None

563 # Numbered objects (1 - 99999)

564 if designation_pf.isdecimal():

565 number = int(designation_pf)

566

567 # Numbered objects (620000+)

568 elif designation_pf[0] == "~":

569 number = 620000

570 number_pf = designation_pf[1:]

571 for i, c in enumerate(number_pf):

572 power = len(number_pf) - (i + 1)

573 number += BASE62_MAP[c] * (62**power)

574

575 # Numbered objects (100000 - 619999)

576 else:

577 number = BASE62_MAP[designation_pf[0]] * 10000 + int(designation_pf[1:])

578

579 if number is None:

580 raise ValueError("Packed numbered designation could not be unpacked.")

581 else:

582 designation = str(number)

583

584 return designation

585

586

587def unpack_provisional_designation(designation_pf: str) -> str:

588 """

589 Unpack a provisional MPC designation.

590

591 Examples of provisional designations:

592 Provisional Unpacked

593 J95X00A 1995 XA

594 J95X01L 1995 XL1

595 J95F13B 1995 FB13

596 J98SA8Q 1998 SQ108

597 J98SC7V 1998 SV127

598 J98SG2S 1998 SS162

599 K99AJ3Z 2099 AZ193

600 K08Aa0A 2008 AA360

601 K07Tf8A 2007 TA418

602

603 Parameters

604 ----------

605 designation_pf : str

606 MPC packed provisional designation.

607

608 Returns

609 -------

610 designation : str

611 MPC unpacked provisional designation.

612

613 Raises

614 ------

615 ValueError : If the provisional designation cannot be unpacked.

616 The packed provisional designation is not 7 characters long.

617 The packed provisional designation does not have a year.

618 """

619 if len(designation_pf) != 7:

620 raise ValueError("Provisional designation must be 7 characters long.")

621 if not designation_pf[1].isdecimal() or not designation_pf[2].isdecimal():

622 raise ValueError("Provisional designation must have a year.")

623 year = str(BASE62_MAP[designation_pf[0]] * 100 + int(designation_pf[1:3]))

624 letter1 = designation_pf[3]

625 letter2 = designation_pf[6]

626 if letter1.isdecimal() or letter2.isdecimal():

627 raise ValueError()

628 cycle1 = designation_pf[4]

629 cycle2 = designation_pf[5]

630

631 number = int(BASE62_MAP[cycle1]) * 10 + BASE62_MAP[cycle2]

632 if number == 0:

633 number_str = ""

634 else:

635 number_str = str(number)

636

637 designation = "{} {}{}{}".format(year, letter1, letter2, number_str)

638

639 return designation

640

641

642def unpack_survey_designation(designation_pf: str) -> str:

643 """

644 Unpack a survey MPC designation.

645

646 Examples of survey designations:

647 Survey Packed

648 PLS2040 2040 P-L

649 T1S3138 3138 T-1

650 T2S1010 1010 T-2

651 T3S4101 4101 T-3

652

653 Parameters

654 ----------

655 designation_pf : str

656 MPC packed survey designation.

657

658 Returns

659 -------

660 designation : str

661 MPC unpacked survey designation.

662

663 Raises

664 ------

665 ValueError : If the survey designation cannot be unpacked.

666 The packed survey designation does not start with PLS, T1S, T2S, or T3S.

667 """

668 number = int(designation_pf[3:8])

669 survey_pf = designation_pf[0:3]

670 if survey_pf not in {"PLS", "T1S", "T2S", "T3S"}:

671 raise ValueError(

672 "Packed survey designation must start with PLS, T1S, T2S, or T3S."

673 )

674

675 if survey_pf == "PLS":

676 survey = "P-L"

677

678 if survey_pf[0] == "T" and survey_pf[2] == "S":

679 survey = "T-{}".format(survey_pf[1])

680

681 designation = "{} {}".format(number, survey)

682 return designation

683

684

685def unpack_mpc_designation(designation_pf: str) -> str:

686 """

687 Unpack a packed MPC designation. For example, provisional

688 designation J98SG2S will be unpacked to 1998 SS162. Permanent

689 designation 00323 will be unpacked to 323.

690

691 TODO: add support for comet and natural satellite designations

692

693 Parameters

694 ----------

695 designation_pf : str

696 MPC packed form designation.

697

698 Returns

699 -------

700 designation : str

701 MPC unpacked designation.

702

703 Raises

704 ------

705 ValueError : If designation_pf cannot be unpacked.

706 """

707 # Lets see if its a numbered object

708 try:

709 return unpack_numbered_designation(designation_pf)

710 except ValueError:

711 pass

712

713 # Lets see if its a provisional designation

714 try:

715 return unpack_provisional_designation(designation_pf)

716 except ValueError:

717 pass

718

719 # Lets see if its a survey designation

720 try:

721 return unpack_survey_designation(designation_pf)

722 except ValueError:

723 pass

724

725 # At this point we haven't had any success so lets raise an error

726 err = (

727 "Packed form designation '{}' could not be unpacked.\n"

728 "It could not be recognized as any of the following:\n"

729 " - a numbered object (e.g. '03202', 'K3289', '~AZaz')\n"

730 " - a provisional designation (e.g. 'J98SC7V', 'K08Aa0A')\n"

731 " - a survey designation (e.g. 'PLS2040', 'T1S3138')"

732 )

733 raise ValueError(err.format(designation_pf))

Coverage for python / lsst / pipe / tasks / associationUtils.py: 9%

207 statements