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:
218 lookup_table_handle = input_ref_dict[
"fgcm_lookup_table"]
222 ref_config.filterMap = self.config.fgcmLoadReferenceCatalog.filterMap
224 for ref
in inputRefs.ref_cat],
225 refCats=butlerQC.get(inputRefs.ref_cat),
226 name=self.config.connections.ref_cat,
229 self.makeSubtask(
'fgcmLoadReferenceCatalog',
230 refObjLoader=ref_obj_loader,
231 refCatName=self.config.connections.ref_cat)
233 lookup_table_handle =
None
239 visit_summary_handle_dict = {handle.dataId[
'visit']: [handle,
None]
for
240 handle
in input_ref_dict[
'visit_summaries']}
242 camera = input_ref_dict[
"camera"]
246 visit_summary_handle_dict=visit_summary_handle_dict,
247 isolated_star_cat_handle_dict=isolated_star_cat_handle_dict,
248 isolated_star_source_handle_dict=isolated_star_source_handle_dict,
249 lookup_table_handle=lookup_table_handle,
252 butlerQC.put(struct.fgcm_visit_catalog, outputRefs.fgcm_visit_catalog)
253 butlerQC.put(struct.fgcm_star_observations, outputRefs.fgcm_star_observations)
254 butlerQC.put(struct.fgcm_star_ids, outputRefs.fgcm_star_ids)
255 if self.config.doReferenceMatches:
256 butlerQC.put(struct.fgcm_reference_stars, outputRefs.fgcm_reference_stars)
258 def run(self, *, camera, visit_summary_handle_dict, isolated_star_cat_handle_dict,
259 isolated_star_source_handle_dict, lookup_table_handle=None):
260 """Run the fgcmBuildFromIsolatedStarsTask.
264 camera : `lsst.afw.cameraGeom.Camera`
266 visit_summary_handle_dict : `dict` [`int`, [`lsst.daf.butler.DeferredDatasetHandle`]]
267 Visit summary dataset handles, with the visit as key.
268 isolated_star_cat_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
269 Isolated star catalog dataset handles, with the tract as key.
270 isolated_star_source_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
271 Isolated star source dataset handles, with the tract as key.
272 lookup_table_handle : `lsst.daf.butler.DeferredDatasetHandle`, optional
273 Data reference to fgcm look-up table (used if matching reference stars).
277 struct : `lsst.pipe.base.struct`
278 Catalogs for persistence, with attributes:
280 ``fgcm_visit_catalog``
281 Catalog of per-visit data (`lsst.afw.table.ExposureCatalog`).
282 ``fgcm_star_observations``
283 Catalog of individual star observations (`astropy.table.Table`).
285 Catalog of per-star ids and positions (`astropy.table.Table`).
286 ``fgcm_reference_stars``
287 Catalog of reference stars matched to fgcm stars (`astropy.table.Table`).
291 calib_flux_aperture_radius =
None
292 if self.config.doSubtractLocalBackground:
294 calib_flux_aperture_radius = computeApertureRadiusFromName(self.config.instFluxField)
295 except RuntimeError
as e:
296 raise RuntimeError(
"Could not determine aperture radius from %s. "
297 "Cannot use doSubtractLocalBackground." %
298 (self.config.instFluxField))
from e
301 if self.config.doReferenceMatches:
302 if lookup_table_handle
is None:
303 raise RuntimeError(
"Must supply lookup_table_handle if config.doReferenceMatches is True.")
309 isolated_star_cat_handle_dict,
310 isolated_star_source_handle_dict,
313 calib_flux_aperture_radius=calib_flux_aperture_radius,
328 if self.config.doReferenceMatches:
333 return pipeBase.Struct(
334 fgcm_visit_catalog=visit_cat,
335 fgcm_star_observations=star_obs,
336 fgcm_star_ids=fgcm_obj,
337 fgcm_reference_stars=fgcm_ref,
342 isolated_star_cat_handle_dict,
343 isolated_star_source_handle_dict,
346 calib_flux_aperture_radius=None,
348 """Make all star observations from isolated star catalogs.
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`
360 calib_flux_aperture_radius : `float`, optional
361 Radius for the calibration flux aperture.
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.
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",
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)
394 if self.sourceSelector.config.doFlags:
395 source_columns.extend(self.sourceSelector.config.flags.bad)
397 if self.config.doSubtractLocalBackground:
398 local_background_area = np.pi*calib_flux_aperture_radius**2.
402 approx_pixel_area_fields = computeApproxPixelAreaFields(camera)
405 detector_mapping = {}
406 for detector_index, detector
in enumerate(camera):
407 detector_mapping[detector.getId()] = detector_index
417 (
"inst_mag_err",
"f4"),
419 (
"delta_mag_bkg",
"f4"),
420 (
"delta_mag_aper",
"f4"),
421 (
"delta_mag_err_aper",
"f4"),
426 (
"isolated_star_id",
"i8"),
429 (
"obs_arr_index",
"i4"),
436 merge_source_counter = 0
440 visit_cat_table = visit_cat.asAstropy()
442 for tract
in sorted(isolated_star_cat_handle_dict):
443 stars = isolated_star_cat_handle_dict[tract].get()
444 sources = isolated_star_source_handle_dict[tract].get(parameters={
"columns": source_columns})
447 good_sources = self.sourceSelector.selectSources(sources).selected
448 if self.config.doSubtractLocalBackground:
449 good_sources &= (~sources[local_background_flag_name])
450 local_background = local_background_area*sources[self.config.localBackgroundFluxField]
451 good_sources &= ((sources[self.config.instFluxField] - local_background) > 0)
453 if good_sources.sum() == 0:
454 self.log.info(
"No good sources found in tract %d", tract)
462 if len(self.config.requiredBands) > 0:
463 loop_bands = self.config.requiredBands
465 loop_bands = np.unique(sources[
"band"])
467 n_req = np.zeros((len(loop_bands), len(stars)), dtype=np.int32)
468 for i, band
in enumerate(loop_bands):
469 (band_use,) = (sources[good_sources][
"band"] == band).nonzero()
472 (i, sources[good_sources][
"obj_index"][band_use]),
476 if len(self.config.requiredBands) > 0:
479 (good_stars,) = (n_req.min(axis=0) >= self.config.minPerBand).nonzero()
483 (good_stars,) = (n_req.max(axis=0) >= self.config.minPerBand).nonzero()
485 if len(good_stars) == 0:
486 self.log.info(
"No good stars found in tract %d", tract)
491 obj_index = sources[
"obj_index"][good_sources]
492 a, b = esutil.numpy_util.match(good_stars, obj_index)
495 _, index_new = np.unique(a, return_index=
True)
496 stars[
"source_cat_index"][good_stars] = index_new
497 sources = sources[good_sources][b]
498 sources[
"obj_index"][:] = a
499 stars = stars[good_stars]
501 nsource = np.zeros(len(stars), dtype=np.int32)
504 sources[
"obj_index"],
507 stars[
"nsource"][:] = nsource
520 star_obs = Table(data=np.zeros(len(sources), dtype=star_obs_dtype))
521 star_obs[
"ra"] = sources[
"ra"]
522 star_obs[
"dec"] = sources[
"dec"]
523 star_obs[
"x"] = sources[
"x"]
524 star_obs[
"y"] = sources[
"y"]
525 star_obs[
"visit"] = sources[
"visit"]
526 star_obs[
"detector"] = sources[
"detector"]
528 visit_match, obs_match = esutil.numpy_util.match(visit_cat_table[
"visit"], sources[
"visit"])
530 exp_time = np.zeros(len(star_obs))
531 exp_time[obs_match] = visit_cat_table[
"exptime"][visit_match]
533 with warnings.catch_warnings():
535 warnings.simplefilter(
"ignore")
537 inst_mag_inner = -2.5*np.log10(sources[self.config.apertureInnerInstFluxField])
538 inst_mag_err_inner = k*(sources[self.config.apertureInnerInstFluxField +
"Err"]
539 / sources[self.config.apertureInnerInstFluxField])
540 inst_mag_outer = -2.5*np.log10(sources[self.config.apertureOuterInstFluxField])
541 inst_mag_err_outer = k*(sources[self.config.apertureOuterInstFluxField +
"Err"]
542 / sources[self.config.apertureOuterInstFluxField])
543 star_obs[
"delta_mag_aper"] = inst_mag_inner - inst_mag_outer
544 star_obs[
"delta_mag_err_aper"] = np.sqrt(inst_mag_err_inner**2. + inst_mag_err_outer**2.)
546 bad = ~np.isfinite(star_obs[
"delta_mag_aper"])
547 star_obs[
"delta_mag_aper"][bad] = 99.0
548 star_obs[
"delta_mag_err_aper"][bad] = 99.0
550 if self.config.doSubtractLocalBackground:
564 local_background = local_background_area*sources[self.config.localBackgroundFluxField]
565 star_obs[
"delta_mag_bkg"] = (-2.5*np.log10(sources[self.config.instFluxField]
566 - local_background) -
567 -2.5*np.log10(sources[self.config.instFluxField]))
570 for detector
in camera:
571 detector_id = detector.getId()
573 (use,) = (star_obs[
"detector"][obs_match] == detector_id).nonzero()
581 jac = approx_pixel_area_fields[detector_id].evaluate(
582 star_obs[
"x"][obs_match][use],
583 star_obs[
"y"][obs_match][use],
585 star_obs[
"jacobian"][obs_match[use]] = jac
586 scaled_inst_flux = (sources[self.config.instFluxField][obs_match[use]]
587 * visit_cat_table[
"scaling"][visit_match[use],
588 detector_mapping[detector_id]])
589 star_obs[
"inst_mag"][obs_match[use]] = (-2.5 * np.log10(scaled_inst_flux
593 star_obs[
"inst_mag_err"] = k*(sources[self.config.instFluxField +
"Err"]
594 / sources[self.config.instFluxField])
597 if self.config.doApplyWcsJacobian:
598 star_obs[
"inst_mag"] -= 2.5*np.log10(star_obs[
"jacobian"])
601 fgcm_obj = Table(data=np.zeros(len(stars), dtype=fgcm_obj_dtype))
602 fgcm_obj[
"isolated_star_id"] = stars[
"isolated_star_id"]
603 fgcm_obj[
"ra"] = stars[
"ra"]
604 fgcm_obj[
"dec"] = stars[
"dec"]
605 fgcm_obj[
"obs_arr_index"] = stars[
"source_cat_index"]
606 fgcm_obj[
"n_obs"] = stars[
"nsource"]
609 fgcm_obj[
"obs_arr_index"] += merge_source_counter
611 fgcm_objs.append(fgcm_obj)
612 star_obs_cats.append(star_obs)
614 merge_source_counter += len(star_obs)
616 fgcm_obj = vstack(fgcm_objs)
619 fgcm_obj[
"fgcm_id"][:] = np.arange(len(fgcm_obj)) + 1
621 return fgcm_obj, vstack(star_obs_cats)
624 """Downsample stars according to density.
626 Catalogs are modified in-place.
630 fgcm_obj : `astropy.table.Table`
631 Catalog of per-star ids and positions.
632 star_obs : `astropy.table.Table`
633 Catalog of individual star observations.
635 if self.config.randomSeed
is not None:
636 np.random.seed(seed=self.config.randomSeed)
638 ipnest = hpg.angle_to_pixel(self.config.densityCutNside, fgcm_obj[
"ra"], fgcm_obj[
"dec"])
643 hist, rev_indices = histogram_rev_sorted(ipnest)
645 obj_use = np.ones(len(fgcm_obj), dtype=bool)
647 (high,) = (hist > self.config.densityCutMaxPerPixel).nonzero()
648 (ok,) = (hist > 0).nonzero()
649 self.log.info(
"There are %d/%d pixels with high stellar density.", high.size, ok.size)
650 for i
in range(high.size):
652 pix_indices = rev_indices[rev_indices[high[i]]: rev_indices[high[i] + 1]]
654 cut = np.random.choice(
656 size=pix_indices.size - self.config.densityCutMaxPerPixel,
661 fgcm_obj = fgcm_obj[obj_use]
663 obs_index = np.zeros(np.sum(fgcm_obj[
"n_obs"]), dtype=np.int32)
665 for i
in range(len(fgcm_obj)):
666 n_obs = fgcm_obj[
"n_obs"][i]
667 obs_index[ctr: ctr + n_obs] = (
668 np.arange(fgcm_obj[
"obs_arr_index"][i], fgcm_obj[
"obs_arr_index"][i] + n_obs)
670 fgcm_obj[
"obs_arr_index"][i] = ctr
673 star_obs = star_obs[obs_index]
690 """Associate fgcm isolated stars with reference stars.
694 lookup_table_handle : `lsst.daf.butler.DeferredDatasetHandle`, optional
695 Data reference to fgcm look-up table (used if matching reference stars).
696 visit_cat : `lsst.afw.table.ExposureCatalog`
697 Catalog of per-visit data.
698 fgcm_obj : `astropy.table.Table`
699 Catalog of per-star ids and positions
703 ref_cat : `astropy.table.Table`
704 Catalog of reference stars matched to fgcm stars.
707 lut_cat = lookup_table_handle.get()
709 std_filter_dict = {filter_name: std_filter
for (filter_name, std_filter)
in
710 zip(lut_cat[0][
"physicalFilters"].split(
","),
711 lut_cat[0][
"stdPhysicalFilters"].split(
","))}
712 std_lambda_dict = {std_filter: std_lambda
for (std_filter, std_lambda)
in
713 zip(lut_cat[0][
"stdPhysicalFilters"].split(
","),
714 lut_cat[0][
"lambdaStdFilter"])}
722 self.log.info(
"Using the following reference filters: %s", (
", ".join(reference_filter_names)))
725 ipnest = hpg.angle_to_pixel(self.config.coarseNside, fgcm_obj[
"ra"], fgcm_obj[
"dec"])
726 hpix, revpix = histogram_rev_sorted(ipnest)
731 dtype = [(
"fgcm_id",
"i4"),
732 (
"refMag",
"f4", (len(reference_filter_names), )),
733 (
"refMagErr",
"f4", (len(reference_filter_names), ))]
735 (gdpix,) = (hpix > 0).nonzero()
736 for ii, gpix
in enumerate(gdpix):
737 p1a = revpix[revpix[gpix]: revpix[gpix + 1]]
740 ra_wrap = fgcm_obj[
"ra"][p1a]
741 if (ra_wrap.min() < 10.0)
and (ra_wrap.max() > 350.0):
742 ra_wrap[ra_wrap > 180.0] -= 360.0
743 mean_ra = np.mean(ra_wrap) % 360.0
745 mean_ra = np.mean(ra_wrap)
746 mean_dec = np.mean(fgcm_obj[
"dec"][p1a])
748 dist = esutil.coords.sphdist(mean_ra, mean_dec, fgcm_obj[
"ra"][p1a], fgcm_obj[
"dec"][p1a])
751 if rad < hpg.nside_to_resolution(self.config.coarseNside)/2.:
757 reference_filter_names,
762 self.config.coarseNside,
763 hpg.nest_to_ring(self.config.coarseNside, ipnest[p1a[0]]),
764 reference_filter_names,
766 if ref_cat.size == 0:
770 with Matcher(fgcm_obj[
"ra"][p1a], fgcm_obj[
"dec"][p1a])
as matcher:
771 idx, i1, i2, d = matcher.query_radius(
774 self.config.matchRadius/3600.,
782 pixel_cat = Table(data=np.zeros(i1.size, dtype=dtype))
783 pixel_cat[
"fgcm_id"] = fgcm_obj[
"fgcm_id"][p1a[i1]]
784 pixel_cat[
"refMag"][:, :] = ref_cat[
"refMag"][i2, :]
785 pixel_cat[
"refMagErr"][:, :] = ref_cat[
"refMagErr"][i2, :]
787 pixel_cats.append(pixel_cat)
790 "Found %d reference matches in pixel %d (%d of %d).",
797 ref_cat_full = vstack(pixel_cats)
798 ref_cat_full.meta.update({
'FILTERNAMES': reference_filter_names})