347 ):
348 """Make all star observations from isolated star catalogs.
349
350 Parameters
351 ----------
352 isolated_star_cat_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
353 Isolated star catalog dataset handles, with the tract as key.
354 isolated_star_source_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
355 Isolated star source dataset handles, with the tract as key.
356 visit_cat : `lsst.afw.table.ExposureCatalog`
357 Catalog of per-visit data.
358 camera : `lsst.afw.cameraGeom.Camera`
359 The camera object.
360 calib_flux_aperture_radius : `float`, optional
361 Radius for the calibration flux aperture.
362
363 Returns
364 -------
365 fgcm_obj : `astropy.table.Table`
366 Catalog of ids and positions for each star.
367 star_obs : `astropy.table.Table`
368 Catalog of individual star observations.
369 """
370 source_columns = [
371 "sourceId",
372 "visit",
373 "detector",
374 "ra",
375 "dec",
376 "x",
377 "y",
378 "physical_filter",
379 "band",
380 "obj_index",
381 self.config.instFluxField,
382 self.config.instFluxField + "Err",
383 self.config.apertureInnerInstFluxField,
384 self.config.apertureInnerInstFluxField + "Err",
385 self.config.apertureOuterInstFluxField,
386 self.config.apertureOuterInstFluxField + "Err",
387 ]
388
389 if self.config.doSubtractLocalBackground:
390 source_columns.append(self.config.localBackgroundFluxField)
391 local_background_flag_name = self.config.localBackgroundFluxField[0: -len('instFlux')] + 'flag'
392 source_columns.append(local_background_flag_name)
393
394 if self.sourceSelector.config.doFlags:
395 source_columns.extend(self.sourceSelector.config.flags.bad)
396
397 local_background_area = np.pi*calib_flux_aperture_radius**2.
398
399
400
401 approx_pixel_area_fields = computeApproxPixelAreaFields(camera)
402
403
404 detector_mapping = {}
405 for detector_index, detector in enumerate(camera):
406 detector_mapping[detector.getId()] = detector_index
407
408 star_obs_dtype = [
409 ("ra", "f8"),
410 ("dec", "f8"),
411 ("x", "f8"),
412 ("y", "f8"),
413 ("visit", "i8"),
414 ("detector", "i4"),
415 ("inst_mag", "f4"),
416 ("inst_mag_err", "f4"),
417 ("jacobian", "f4"),
418 ("delta_mag_bkg", "f4"),
419 ("delta_mag_aper", "f4"),
420 ("delta_mag_err_aper", "f4"),
421 ]
422
423 fgcm_obj_dtype = [
424 ("fgcm_id", "i8"),
425 ("isolated_star_id", "i8"),
426 ("ra", "f8"),
427 ("dec", "f8"),
428 ("obs_arr_index", "i4"),
429 ("n_obs", "i4"),
430 ("obj_flag", "i4"),
431 ]
432
433 fgcm_objs = []
434 star_obs_cats = []
435 merge_source_counter = 0
436
437 k = 2.5/np.log(10.)
438
439 visit_cat_table = visit_cat.asAstropy()
440
441 for tract in sorted(isolated_star_cat_handle_dict):
442 stars = isolated_star_cat_handle_dict[tract].get()
443 sources = isolated_star_source_handle_dict[tract].get(parameters={"columns": source_columns})
444
445
446 good_sources = self.sourceSelector.selectSources(sources).selected
447 if self.config.doSubtractLocalBackground:
448 good_sources &= (~sources[local_background_flag_name])
449 local_background = local_background_area*sources[self.config.localBackgroundFluxField]
450 good_sources &= ((sources[self.config.instFluxField] - local_background) > 0)
451
452 if good_sources.sum() == 0:
453 self.log.info("No good sources found in tract %d", tract)
454 continue
455
456
457
458
459
460
461 if len(self.config.requiredBands) > 0:
462 loop_bands = self.config.requiredBands
463 else:
464 loop_bands = np.unique(sources["band"])
465
466 n_req = np.zeros((len(loop_bands), len(stars)), dtype=np.int32)
467 for i, band in enumerate(loop_bands):
468 (band_use,) = (sources[good_sources]["band"] == band).nonzero()
469 np.add.at(
470 n_req,
471 (i, sources[good_sources]["obj_index"][band_use]),
472 1,
473 )
474
475 if len(self.config.requiredBands) > 0:
476
477
478 (good_stars,) = (n_req.min(axis=0) >= self.config.minPerBand).nonzero()
479 else:
480
481
482 (good_stars,) = (n_req.max(axis=0) >= self.config.minPerBand).nonzero()
483
484 if len(good_stars) == 0:
485 self.log.info("No good stars found in tract %d", tract)
486 continue
487
488
489
490 obj_index = sources["obj_index"][good_sources]
491 a, b = esutil.numpy_util.match(good_stars, obj_index)
492
493
494 _, index_new = np.unique(a, return_index=True)
495 stars["source_cat_index"][good_stars] = index_new
496 sources = sources[good_sources][b]
497 sources["obj_index"][:] = a
498 stars = stars[good_stars]
499
500 nsource = np.zeros(len(stars), dtype=np.int32)
501 np.add.at(
502 nsource,
503 sources["obj_index"],
504 1,
505 )
506 stars["nsource"][:] = nsource
507
508
509
510
511
512
513
514
515
516
517
518
519 star_obs = Table(data=np.zeros(len(sources), dtype=star_obs_dtype))
520 star_obs["ra"] = sources["ra"]
521 star_obs["dec"] = sources["dec"]
522 star_obs["x"] = sources["x"]
523 star_obs["y"] = sources["y"]
524 star_obs["visit"] = sources["visit"]
525 star_obs["detector"] = sources["detector"]
526
527 visit_match, obs_match = esutil.numpy_util.match(visit_cat_table["visit"], sources["visit"])
528
529 exp_time = np.zeros(len(star_obs))
530 exp_time[obs_match] = visit_cat_table["exptime"][visit_match]
531
532 with warnings.catch_warnings():
533
534 warnings.simplefilter("ignore")
535
536 inst_mag_inner = -2.5*np.log10(sources[self.config.apertureInnerInstFluxField])
537 inst_mag_err_inner = k*(sources[self.config.apertureInnerInstFluxField + "Err"]
538 / sources[self.config.apertureInnerInstFluxField])
539 inst_mag_outer = -2.5*np.log10(sources[self.config.apertureOuterInstFluxField])
540 inst_mag_err_outer = k*(sources[self.config.apertureOuterInstFluxField + "Err"]
541 / sources[self.config.apertureOuterInstFluxField])
542 star_obs["delta_mag_aper"] = inst_mag_inner - inst_mag_outer
543 star_obs["delta_mag_err_aper"] = np.sqrt(inst_mag_err_inner**2. + inst_mag_err_outer**2.)
544
545 bad = ~np.isfinite(star_obs["delta_mag_aper"])
546 star_obs["delta_mag_aper"][bad] = 99.0
547 star_obs["delta_mag_err_aper"][bad] = 99.0
548
549 if self.config.doSubtractLocalBackground:
550
551
552
553
554
555
556
557
558
559
560
561
562
563 local_background = local_background_area*sources[self.config.localBackgroundFluxField]
564 star_obs["delta_mag_bkg"] = (-2.5*np.log10(sources[self.config.instFluxField]
565 - local_background) -
566 -2.5*np.log10(sources[self.config.instFluxField]))
567
568
569 for detector in camera:
570 detector_id = detector.getId()
571
572 (use,) = (star_obs["detector"][obs_match] == detector_id).nonzero()
573
574
575
576
577
578
579
580 jac = approx_pixel_area_fields[detector_id].evaluate(
581 star_obs["x"][obs_match][use],
582 star_obs["y"][obs_match][use],
583 )
584 star_obs["jacobian"][obs_match[use]] = jac
585 scaled_inst_flux = (sources[self.config.instFluxField][obs_match[use]]
586 * visit_cat_table["scaling"][visit_match[use],
587 detector_mapping[detector_id]])
588 star_obs["inst_mag"][obs_match[use]] = (-2.5 * np.log10(scaled_inst_flux
589 / exp_time[use]))
590
591
592 star_obs["inst_mag_err"] = k*(sources[self.config.instFluxField + "Err"]
593 / sources[self.config.instFluxField])
594
595
596 if self.config.doApplyWcsJacobian:
597 star_obs["inst_mag"] -= 2.5*np.log10(star_obs["jacobian"])
598
599
600 fgcm_obj = Table(data=np.zeros(len(stars), dtype=fgcm_obj_dtype))
601 fgcm_obj["isolated_star_id"] = stars["isolated_star_id"]
602 fgcm_obj["ra"] = stars["ra"]
603 fgcm_obj["dec"] = stars["dec"]
604 fgcm_obj["obs_arr_index"] = stars["source_cat_index"]
605 fgcm_obj["n_obs"] = stars["nsource"]
606
607
608 fgcm_obj["obs_arr_index"] += merge_source_counter
609
610 fgcm_objs.append(fgcm_obj)
611 star_obs_cats.append(star_obs)
612
613 merge_source_counter += len(star_obs)
614
615 fgcm_obj = vstack(fgcm_objs)
616
617
618 fgcm_obj["fgcm_id"][:] = np.arange(len(fgcm_obj)) + 1
619
620 return fgcm_obj, vstack(star_obs_cats)
621