198 input_ref_dict = butlerQC.get(inputRefs)
200 isolated_star_cat_handles = input_ref_dict[
"isolated_star_cats"]
201 isolated_star_source_handles = input_ref_dict[
"isolated_star_sources"]
203 isolated_star_cat_handle_dict = {
204 handle.dataId[
"tract"]: handle
for handle
in isolated_star_cat_handles
206 isolated_star_source_handle_dict = {
207 handle.dataId[
"tract"]: handle
for handle
in isolated_star_source_handles
210 if len(isolated_star_cat_handle_dict) != len(isolated_star_source_handle_dict):
211 raise RuntimeError(
"isolated_star_cats and isolate_star_sources must have same length.")
213 for tract
in isolated_star_cat_handle_dict:
214 if tract
not in isolated_star_source_handle_dict:
215 raise RuntimeError(f
"tract {tract} in isolated_star_cats but not isolated_star_sources")
217 if self.config.doReferenceMatches:
221 ref_config.filterMap = self.config.fgcmLoadReferenceCatalog.filterMap
223 for ref
in inputRefs.ref_cat],
224 refCats=butlerQC.get(inputRefs.ref_cat),
225 name=self.config.connections.ref_cat,
228 self.makeSubtask(
'fgcmLoadReferenceCatalog',
229 refObjLoader=ref_obj_loader,
230 refCatName=self.config.connections.ref_cat)
236 visit_summary_handle_dict = {handle.dataId[
'visit']: [handle,
None]
for
237 handle
in input_ref_dict[
'visit_summaries']}
239 camera = input_ref_dict[
"camera"]
241 lookup_table_handle = input_ref_dict[
"fgcm_lookup_table"]
245 visit_summary_handle_dict=visit_summary_handle_dict,
246 isolated_star_cat_handle_dict=isolated_star_cat_handle_dict,
247 isolated_star_source_handle_dict=isolated_star_source_handle_dict,
248 lookup_table_handle=lookup_table_handle,
251 butlerQC.put(struct.fgcm_visit_catalog, outputRefs.fgcm_visit_catalog)
252 butlerQC.put(struct.fgcm_star_observations, outputRefs.fgcm_star_observations)
253 butlerQC.put(struct.fgcm_star_ids, outputRefs.fgcm_star_ids)
254 if self.config.doReferenceMatches:
255 butlerQC.put(struct.fgcm_reference_stars, outputRefs.fgcm_reference_stars)
257 def run(self, *, camera, visit_summary_handle_dict, isolated_star_cat_handle_dict,
258 isolated_star_source_handle_dict, lookup_table_handle):
259 """Run the fgcmBuildFromIsolatedStarsTask.
263 camera : `lsst.afw.cameraGeom.Camera`
265 visit_summary_handle_dict : `dict` [`int`, [`lsst.daf.butler.DeferredDatasetHandle`]]
266 Visit summary dataset handles, with the visit as key.
267 isolated_star_cat_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
268 Isolated star catalog dataset handles, with the tract as key.
269 isolated_star_source_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
270 Isolated star source dataset handles, with the tract as key.
271 lookup_table_handle : `lsst.daf.butler.DeferredDatasetHandle`
272 Data reference to fgcm look-up table.
276 struct : `lsst.pipe.base.struct`
277 Catalogs for persistence, with attributes:
279 ``fgcm_visit_catalog``
280 Catalog of per-visit data (`lsst.afw.table.ExposureCatalog`).
281 ``fgcm_star_observations``
282 Catalog of individual star observations (`astropy.table.Table`).
284 Catalog of per-star ids and positions (`astropy.table.Table`).
285 ``fgcm_reference_stars``
286 Catalog of reference stars matched to fgcm stars (`astropy.table.Table`).
290 calib_flux_aperture_radius =
None
291 if self.config.doSubtractLocalBackground:
293 calib_flux_aperture_radius = computeApertureRadiusFromName(self.config.instFluxField)
294 except RuntimeError
as e:
295 raise RuntimeError(
"Could not determine aperture radius from %s. "
296 "Cannot use doSubtractLocalBackground." %
297 (self.config.instFluxField))
from e
300 if self.config.doReferenceMatches:
301 if lookup_table_handle
is None:
302 raise RuntimeError(
"Must supply lookup_table_handle if config.doReferenceMatches is True.")
304 lut_cat = lookup_table_handle.get()
305 if len(camera) == lut_cat[0][
"nCcd"]:
306 use_science_detectors =
False
313 use_science_detectors =
True
318 visit_summary_handle_dict,
319 useScienceDetectors=use_science_detectors,
324 isolated_star_cat_handle_dict,
325 isolated_star_source_handle_dict,
328 calib_flux_aperture_radius=calib_flux_aperture_radius,
329 use_science_detectors=use_science_detectors,
344 if self.config.doReferenceMatches:
349 return pipeBase.Struct(
350 fgcm_visit_catalog=visit_cat,
351 fgcm_star_observations=star_obs,
352 fgcm_star_ids=fgcm_obj,
353 fgcm_reference_stars=fgcm_ref,
358 isolated_star_cat_handle_dict,
359 isolated_star_source_handle_dict,
362 calib_flux_aperture_radius=None,
363 use_science_detectors=False,
365 """Make all star observations from isolated star catalogs.
369 isolated_star_cat_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
370 Isolated star catalog dataset handles, with the tract as key.
371 isolated_star_source_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
372 Isolated star source dataset handles, with the tract as key.
373 visit_cat : `lsst.afw.table.ExposureCatalog`
374 Catalog of per-visit data.
375 camera : `lsst.afw.cameraGeom.Camera`
377 calib_flux_aperture_radius : `float`, optional
378 Radius for the calibration flux aperture.
379 use_science_detectors : `bool`, optional
380 Only use science detectors?
384 fgcm_obj : `astropy.table.Table`
385 Catalog of ids and positions for each star.
386 star_obs : `astropy.table.Table`
387 Catalog of individual star observations.
400 self.config.instFluxField,
401 self.config.instFluxField +
"Err",
402 self.config.apertureInnerInstFluxField,
403 self.config.apertureInnerInstFluxField +
"Err",
404 self.config.apertureOuterInstFluxField,
405 self.config.apertureOuterInstFluxField +
"Err",
408 if self.config.doSubtractLocalBackground:
409 source_columns.append(self.config.localBackgroundFluxField)
410 local_background_flag_name = self.config.localBackgroundFluxField[0: -len(
'instFlux')] +
'flag'
411 source_columns.append(local_background_flag_name)
413 if self.sourceSelector.config.doFlags:
414 source_columns.extend(self.sourceSelector.config.flags.bad)
416 if self.config.doSubtractLocalBackground:
417 local_background_area = np.pi*calib_flux_aperture_radius**2.
421 approx_pixel_area_fields = computeApproxPixelAreaFields(camera)
424 detector_mapping = {}
425 for detector_index, detector
in enumerate(camera):
426 detector_mapping[detector.getId()] = detector_index
436 (
"inst_mag_err",
"f4"),
438 (
"delta_mag_bkg",
"f4"),
439 (
"delta_mag_aper",
"f4"),
440 (
"delta_mag_err_aper",
"f4"),
445 (
"isolated_star_id",
"i8"),
448 (
"obs_arr_index",
"i4"),
455 merge_source_counter = 0
459 visit_cat_table = visit_cat.asAstropy()
461 for tract
in sorted(isolated_star_cat_handle_dict):
462 stars = isolated_star_cat_handle_dict[tract].get()
463 sources = isolated_star_source_handle_dict[tract].get(parameters={
"columns": source_columns})
466 good_sources = self.sourceSelector.selectSources(sources).selected
467 if self.config.doSubtractLocalBackground:
468 good_sources &= (~sources[local_background_flag_name])
469 local_background = local_background_area*sources[self.config.localBackgroundFluxField]
470 good_sources &= ((sources[self.config.instFluxField] - local_background) > 0)
474 if good_sources.sum() == 0:
475 self.log.info(
"No good sources found in tract %d", tract)
482 for detector
in camera:
483 if use_science_detectors:
484 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE:
486 detector_ids.append(detector.getId())
488 missing_sources1 = np.ones(len(sources), dtype=np.bool_)
489 _, sources_match = esutil.numpy_util.match(visit_cat_table[
"visit"], sources[
"visit"])
490 missing_sources1[sources_match] =
False
491 good_sources &= (~missing_sources1)
493 missing_sources2 = np.ones(len(sources), dtype=np.bool_)
494 _, sources_match = esutil.numpy_util.match(detector_ids, sources[
"detector"])
495 missing_sources2[sources_match] =
False
496 good_sources &= (~missing_sources2)
499 if good_sources.sum() == 0:
500 self.log.info(
"No good sources found in tract %d", tract)
503 self.log.info(
"Tract %d contains %d good sources.", tract, good_sources.sum())
510 if len(self.config.requiredBands) > 0:
511 loop_bands = self.config.requiredBands
513 loop_bands = np.unique(sources[
"band"])
515 n_req = np.zeros((len(loop_bands), len(stars)), dtype=np.int32)
516 for i, band
in enumerate(loop_bands):
517 (band_use,) = (sources[good_sources][
"band"] == band).nonzero()
520 (i, sources[good_sources][
"obj_index"][band_use]),
524 if len(self.config.requiredBands) > 0:
527 (good_stars,) = (n_req.min(axis=0) >= self.config.minPerBand).nonzero()
531 (good_stars,) = (n_req.max(axis=0) >= self.config.minPerBand).nonzero()
533 if len(good_stars) == 0:
534 self.log.info(
"No good stars found in tract %d", tract)
539 obj_index = sources[
"obj_index"][good_sources]
540 a, b = esutil.numpy_util.match(good_stars, obj_index)
543 _, index_new = np.unique(a, return_index=
True)
544 stars[
"source_cat_index"][good_stars] = index_new
545 sources = sources[good_sources][b]
546 sources[
"obj_index"][:] = a
547 stars = stars[good_stars]
549 nsource = np.zeros(len(stars), dtype=np.int32)
552 sources[
"obj_index"],
555 stars[
"nsource"][:] = nsource
568 star_obs = Table(data=np.zeros(len(sources), dtype=star_obs_dtype))
569 star_obs[
"ra"] = sources[
"ra"]
570 star_obs[
"dec"] = sources[
"dec"]
571 star_obs[
"x"] = sources[
"x"]
572 star_obs[
"y"] = sources[
"y"]
573 star_obs[
"visit"] = sources[
"visit"]
574 star_obs[
"detector"] = sources[
"detector"]
576 visit_match, obs_match = esutil.numpy_util.match(visit_cat_table[
"visit"], sources[
"visit"])
578 exp_time = np.zeros(len(star_obs))
579 exp_time[obs_match] = visit_cat_table[
"exptime"][visit_match]
581 with warnings.catch_warnings():
583 warnings.simplefilter(
"ignore")
585 inst_mag_inner = -2.5*np.log10(sources[self.config.apertureInnerInstFluxField])
586 inst_mag_err_inner = k*(sources[self.config.apertureInnerInstFluxField +
"Err"]
587 / sources[self.config.apertureInnerInstFluxField])
588 inst_mag_outer = -2.5*np.log10(sources[self.config.apertureOuterInstFluxField])
589 inst_mag_err_outer = k*(sources[self.config.apertureOuterInstFluxField +
"Err"]
590 / sources[self.config.apertureOuterInstFluxField])
591 star_obs[
"delta_mag_aper"] = inst_mag_inner - inst_mag_outer
592 star_obs[
"delta_mag_err_aper"] = np.sqrt(inst_mag_err_inner**2. + inst_mag_err_outer**2.)
594 bad = ~np.isfinite(star_obs[
"delta_mag_aper"])
595 star_obs[
"delta_mag_aper"][bad] = 99.0
596 star_obs[
"delta_mag_err_aper"][bad] = 99.0
598 if self.config.doSubtractLocalBackground:
612 local_background = local_background_area*sources[self.config.localBackgroundFluxField]
613 star_obs[
"delta_mag_bkg"] = (-2.5*np.log10(sources[self.config.instFluxField]
614 - local_background) -
615 -2.5*np.log10(sources[self.config.instFluxField]))
618 for detector
in camera:
619 if use_science_detectors:
620 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE:
623 detector_id = detector.getId()
625 (use,) = (star_obs[
"detector"][obs_match] == detector_id).nonzero()
633 jac = approx_pixel_area_fields[detector_id].evaluate(
634 star_obs[
"x"][obs_match][use],
635 star_obs[
"y"][obs_match][use],
637 star_obs[
"jacobian"][obs_match[use]] = jac
638 scaled_inst_flux = (sources[self.config.instFluxField][obs_match[use]]
639 * visit_cat_table[
"scaling"][visit_match[use],
640 detector_mapping[detector_id]])
641 star_obs[
"inst_mag"][obs_match[use]] = (-2.5 * np.log10(scaled_inst_flux
647 star_obs[
"inst_mag_err"] = k*(sources[self.config.instFluxField +
"Err"]
648 / sources[self.config.instFluxField])
651 if self.config.doApplyWcsJacobian:
652 star_obs[
"inst_mag"] -= 2.5*np.log10(star_obs[
"jacobian"])
655 fgcm_obj = Table(data=np.zeros(len(stars), dtype=fgcm_obj_dtype))
656 fgcm_obj[
"isolated_star_id"] = stars[
"isolated_star_id"]
657 fgcm_obj[
"ra"] = stars[
"ra"]
658 fgcm_obj[
"dec"] = stars[
"dec"]
659 fgcm_obj[
"obs_arr_index"] = stars[
"source_cat_index"]
660 fgcm_obj[
"n_obs"] = stars[
"nsource"]
663 fgcm_obj[
"obs_arr_index"] += merge_source_counter
665 fgcm_objs.append(fgcm_obj)
666 star_obs_cats.append(star_obs)
668 merge_source_counter += len(star_obs)
670 fgcm_obj = vstack(fgcm_objs)
673 fgcm_obj[
"fgcm_id"][:] = np.arange(len(fgcm_obj)) + 1
675 return fgcm_obj, vstack(star_obs_cats)
744 """Associate fgcm isolated stars with reference stars.
748 lookup_table_handle : `lsst.daf.butler.DeferredDatasetHandle`, optional
749 Data reference to fgcm look-up table (used if matching reference stars).
750 visit_cat : `lsst.afw.table.ExposureCatalog`
751 Catalog of per-visit data.
752 fgcm_obj : `astropy.table.Table`
753 Catalog of per-star ids and positions
757 ref_cat : `astropy.table.Table`
758 Catalog of reference stars matched to fgcm stars.
761 lut_cat = lookup_table_handle.get()
763 std_filter_dict = {filter_name: std_filter
for (filter_name, std_filter)
in
764 zip(lut_cat[0][
"physicalFilters"].split(
","),
765 lut_cat[0][
"stdPhysicalFilters"].split(
","))}
766 std_lambda_dict = {std_filter: std_lambda
for (std_filter, std_lambda)
in
767 zip(lut_cat[0][
"stdPhysicalFilters"].split(
","),
768 lut_cat[0][
"lambdaStdFilter"])}
776 self.log.info(
"Using the following reference filters: %s", (
", ".join(reference_filter_names)))
779 ipnest = hpg.angle_to_pixel(self.config.coarseNside, fgcm_obj[
"ra"], fgcm_obj[
"dec"])
780 hpix, revpix = histogram_rev_sorted(ipnest)
785 dtype = [(
"fgcm_id",
"i4"),
786 (
"refMag",
"f4", (len(reference_filter_names), )),
787 (
"refMagErr",
"f4", (len(reference_filter_names), ))]
789 (gdpix,) = (hpix > 0).nonzero()
790 for ii, gpix
in enumerate(gdpix):
791 p1a = revpix[revpix[gpix]: revpix[gpix + 1]]
794 ra_wrap = fgcm_obj[
"ra"][p1a]
795 if (ra_wrap.min() < 10.0)
and (ra_wrap.max() > 350.0):
796 ra_wrap[ra_wrap > 180.0] -= 360.0
797 mean_ra = np.mean(ra_wrap) % 360.0
799 mean_ra = np.mean(ra_wrap)
800 mean_dec = np.mean(fgcm_obj[
"dec"][p1a])
802 dist = esutil.coords.sphdist(mean_ra, mean_dec, fgcm_obj[
"ra"][p1a], fgcm_obj[
"dec"][p1a])
805 if rad < hpg.nside_to_resolution(self.config.coarseNside)/2.:
811 reference_filter_names,
816 self.config.coarseNside,
817 hpg.nest_to_ring(self.config.coarseNside, ipnest[p1a[0]]),
818 reference_filter_names,
820 if ref_cat.size == 0:
824 with Matcher(fgcm_obj[
"ra"][p1a], fgcm_obj[
"dec"][p1a])
as matcher:
825 idx, i1, i2, d = matcher.query_radius(
828 self.config.matchRadius/3600.,
836 pixel_cat = Table(data=np.zeros(i1.size, dtype=dtype))
837 pixel_cat[
"fgcm_id"] = fgcm_obj[
"fgcm_id"][p1a[i1]]
838 pixel_cat[
"refMag"][:, :] = ref_cat[
"refMag"][i2, :]
839 pixel_cat[
"refMagErr"][:, :] = ref_cat[
"refMagErr"][i2, :]
841 pixel_cats.append(pixel_cat)
844 "Found %d reference matches in pixel %d (%d of %d).",
851 ref_cat_full = vstack(pixel_cats)
852 ref_cat_full.meta.update({
'FILTERNAMES': reference_filter_names})
857 """Compute delta aperture summary statistics (per visit).
861 visit_cat : `lsst.afw.table.ExposureCatalog`
862 Catalog of per-visit data.
863 star_obs : `astropy.table.Table`
864 Catalog of individual star observations.
866 (ok,) = ((star_obs[
"delta_mag_aper"] < 99.0)
867 & (star_obs[
"delta_mag_err_aper"] < 99.0)).nonzero()
869 visit_index = np.zeros(len(star_obs[ok]), dtype=np.int32)
870 a, b = esutil.numpy_util.match(visit_cat[
"visit"], star_obs[
"visit"][ok])
873 h, rev = histogram_rev_sorted(visit_index)
875 visit_cat[
"deltaAper"] = -9999.0
876 h_use, = np.where(h >= 3)
878 i1a = rev[rev[index]: rev[index + 1]]
879 visit_cat[
"deltaAper"][index] = np.median(np.asarray(star_obs[
"delta_mag_aper"][ok[i1a]]))
881 n_detector = visit_cat.schema[
"deltaAperDetector"].asKey().getSize()
885 visit_detector_hash = visit_index * (n_detector + 1) + star_obs[
"detector"][ok]
887 h, rev = histogram_rev_sorted(visit_detector_hash)
889 visit_cat[
"deltaAperDetector"][:] = -9999.0
890 h_use, = np.where(h >= 3)
892 i1a = rev[rev[index]: rev[index + 1]]
893 v_ind = visit_index[i1a[0]]
894 d_ind = visit_detector_hash[i1a[0]] % (n_detector + 1)
895 delta_mag = np.median(np.asarray(star_obs[
"delta_mag_aper"][ok[i1a]]))
896 visit_cat[
"deltaAperDetector"][v_ind, d_ind] = delta_mag