23 Statistics of jointcal vs. single-frame procesing and diagnostic plots.
27 Some of the algorithms and data structures in this code are temporary
28 kludges and will no longer be necessary once the following are available:
30 - a composite data structure that contains all ccds from a single visit
31 - an n-way matching system that preserves the separations between sources
37 from astropy
import units
as u
44 __all__ = [
'JointcalStatistics']
46 MatchDict = collections.namedtuple(
'MatchDict', [
'relative',
'absolute'])
52 Compute statistics on jointcal-processed data, and optionally generate plots.
56 Instantiate JointcalStatistics and call compute_rms() to get the relevant
57 statistics for e.g. unittests, and call make_plots() to generate a suite of
61 def __init__(self, match_radius=0.1*arcseconds, flux_limit=100.0,
62 do_photometry=True, do_astrometry=True,
67 match_radius : lsst.geom.Angle
68 match sources within this radius for RMS statistics
70 Signal/Noise (flux/fluxErr) for sources to be included in the RMS cross-match.
71 100 is a balance between good centroids and enough sources.
72 do_photometry : bool, optional
73 Perform calculations/make plots for photometric metrics.
74 do_astrometry : bool, optional
75 Perform calculations/make plots for astrometric metrics.
76 verbose : bool, optional
84 self.
loglog = _LOG.getChild(
'JointcalStatistics')
88 Match all data_refs to compute the RMS, for all detections above self.flux_limit.
92 data_refs : list of lsst.daf.persistence.butlerSubset.ButlerDataRef
93 A list of data refs to do the calculations between.
94 reference : lsst reference catalog
95 reference catalog to do absolute matching against.
101 Post-jointcal relative RMS of the matched sources.
103 Post-jointcal absolute RMS of matched sources.
105 Post-jointcal photometric repeatability (PA1 from the SRD).
109 self.
filtersfilters = [ref.get(
'calexp_filterLabel').bandLabel
for ref
in data_refs]
112 def compute(catalogs, photoCalibs):
113 """Compute the relative and absolute matches in distance and flux."""
123 refcalib = photoCalibs[0]
if photoCalibs != []
else None
124 dist_rel, flux_rel, ref_flux_rel, source_rel = self.
_make_match_dict_make_match_dict(refcat,
128 dist_abs, flux_abs, ref_flux_abs, source_abs = self.
_make_match_dict_make_match_dict(reference,
133 ref_flux =
MatchDict(ref_flux_rel, ref_flux_abs)
134 source =
MatchDict(source_rel, source_abs)
135 return dist, flux, ref_flux, source
137 old_cats = [ref.get(
'src')
for ref
in data_refs]
142 old_calibs = [ref.get(
'calexp_photoCalib')
for ref
in data_refs]
144 self.old_dist, self.old_flux, self.old_ref_flux, self.
old_sourceold_source = compute(old_cats, old_calibs)
147 new_cats = [ref.get(
'src')
for ref
in data_refs]
150 new_wcss = [ref.get(
'jointcal_wcs')
for ref
in data_refs]
153 new_calibs = [ref.get(
'jointcal_photoCalib')
for ref
in data_refs]
155 for wcs, cat
in zip(new_wcss, new_cats):
159 self.new_dist, self.new_flux, self.new_ref_flux, self.
new_sourcenew_source = compute(new_cats, new_calibs)
162 print(
'old, new relative distance matches:',
163 len(self.old_dist.relative), len(self.new_dist.relative))
164 print(
'old, new absolute distance matches:',
165 len(self.old_dist.absolute), len(self.new_dist.absolute))
166 print(
'old, new relative flux matches:',
167 len(self.old_flux.relative), len(self.new_flux.relative))
168 print(
'old, new absolute flux matches:',
169 len(self.old_flux.absolute), len(self.new_flux.absolute))
177 """Compute the total rms across all sources."""
178 total = sum(sum(dd**2)
for dd
in data.values())
179 n = sum(len(dd)
for dd
in data.values())
180 return np.sqrt(total/n)
191 print(
"relative, absolute astrometry:",
195 print(
"Photometric Accuracy (PA1):", self.
new_PA1new_PA1)
196 Rms_result = collections.namedtuple(
"Rms_result", [
"dist_relative",
"dist_absolute",
"pa1"])
200 name='', interactive=False, per_ccd_plot=False, outdir='.plots'):
202 Make plots of various quantites to help with debugging.
203 Requires that `compute_rms()` was run first.
207 data_refs : list of lsst.daf.persistence.butlerSubset.ButlerDataRef
208 A list of data refs to do the calculations between.
209 old_wcs_list : list of lsst.afw.image.wcs.Wcs
210 A list of the old (pre-jointcal) WCSs, one-to-one corresponding to data_refs.
212 Name to include in plot titles and save files.
214 Turn on matplotlib interactive mode and drop into a debugger when
215 plotting is finished. Otherwise, use a non-interactive backend.
217 Plot the WCS per CCD (takes longer and generates many plots for a large camera)
219 directory to save plots to.
225 matplotlib.use(
'pdf')
227 import matplotlib.pyplot
as plt
228 import astropy.visualization
230 astropy.visualization.quantity_support()
234 self.
loglog.
info(
"N data_refs: %d", len(data_refs))
240 name=name, outdir=outdir)
244 def rms_per_source(data):
245 """Each element of data must already be the "delta" of whatever measurement."""
246 return (np.sqrt([np.mean(dd**2)
for dd
in data.values()])*u.radian).
to(u.arcsecond)
257 new_dist_rms.relative, new_dist_rms.absolute,
260 name=name, outdir=outdir)
263 per_ccd_plot=per_ccd_plot,
264 name=name, outdir=outdir)
271 def _photometric_rms(self, sn_cut=300, magnitude_range=3):
273 Compute the photometric RMS and the photometric repeatablity values (PA1).
278 The minimum signal/noise for sources to be included in the PA1 calculation.
279 magnitude_range : float
280 The range of magnitudes above sn_cut to include in the PA1 calculation.
282 def rms(flux, ref_flux):
283 return np.sqrt([np.mean((ref_flux[dd] - flux[dd])**2)
for dd
in flux])
289 self.
old_refold_ref = np.fromiter(self.old_ref_flux.absolute.values(), dtype=float)
290 self.
new_refnew_ref = np.fromiter(self.new_ref_flux.absolute.values(), dtype=float)
294 def signal_to_noise(sources, flux_key='slot_PsfFlux_instFlux', sigma_key='slot_PsfFlux_instFluxErr'):
295 """Compute the mean signal/noise per source from a MatchDict of SourceRecords."""
296 result = np.empty(len(sources))
297 for i, src
in enumerate(sources.values()):
298 result[i] = np.mean([x[flux_key]/x[sigma_key]
for x
in src])
301 old_sn = signal_to_noise(self.
old_sourceold_source.absolute)
306 print(
"PA1 Magnitude range: {:.3f}, {:.3f}".
format(self.
brightbright, self.
faintfaint))
314 def _make_match_dict(self, reference, visit_catalogs, photoCalibs, refcalib=None):
316 Return several dicts of sourceID:[values] over the catalogs, to be used in RMS calculations.
320 reference : lsst.afw.table.SourceCatalog
321 Catalog to do the matching against.
322 visit_catalogs : list of lsst.afw.table.SourceCatalog
323 Visit source catalogs (values() produced by _make_visit_catalogs)
324 to cross-match against reference.
325 photoCalibs : list of lsst.afw.image.PhotoCalib
326 Exposure PhotoCalibs, 1-1 coorespondent with visit_catalogs.
327 refcalib : lsst.afw.image.PhotoCalib or None
328 Pass a PhotoCalib here to use it to compute nanojansky from the
329 reference catalog ADU slot_flux.
334 dict of sourceID: array(separation distances for that source)
336 dict of sourceID: array(fluxes (nJy) for that source)
338 dict of sourceID: flux (nJy) of the reference object
340 dict of sourceID: list(each SourceRecord that was position-matched
344 if photoCalibs == []:
345 photoCalibs = [[]]*len(visit_catalogs)
347 distances = collections.defaultdict(list)
348 fluxes = collections.defaultdict(list)
350 sources = collections.defaultdict(list)
351 if 'slot_CalibFlux_instFlux' in reference.schema:
352 ref_flux_key =
'slot_CalibFlux'
354 ref_flux_key =
'{}_flux'
356 def get_fluxes(photoCalib, match):
357 """Return (flux, ref_flux) or None if either is invalid."""
359 flux = match[1][
'slot_CalibFlux_instFlux']
363 flux = photoCalib.instFluxToNanojansky(match[1],
"slot_CalibFlux").value
366 if 'slot' in ref_flux_key:
367 ref_flux = match[0][ref_flux_key+
'_instFlux']
371 ref_flux = refcalib.instFluxToNanojansky(match[0], ref_flux_key).value
374 ref_flux = match[0][ref_flux_key.format(filt)]
378 Flux = collections.namedtuple(
'Flux', (
'flux',
'ref_flux'))
379 return Flux(flux, ref_flux)
381 for cat, photoCalib, filt
in zip(visit_catalogs, photoCalibs, self.
filtersfilters):
382 good = (cat.get(
'base_PsfFlux_instFlux')/cat.get(
'base_PsfFlux_instFluxErr')) > self.
flux_limitflux_limit
384 good &= (cat.get(
'base_ClassificationExtendedness_value') == 0)
388 flux = get_fluxes(photoCalib, m)
392 fluxes[m[0].getId()].
append(flux.flux)
394 ref_fluxes[m[0].getId()] = flux.ref_flux
398 distances[m[0].getId()].
append(m[2])
400 sources[m[0].getId()].
append(m[1])
402 for source
in distances:
403 distances[source] = np.array(distances[source])
404 for source
in fluxes:
405 fluxes[source] = np.array(fluxes[source])
407 return distances, fluxes, ref_fluxes, sources
409 def _make_visit_catalogs(self, catalogs, visits):
411 Merge all catalogs from the each visit.
412 NOTE: creating this structure is somewhat slow, and will be unnecessary
413 once a full-visit composite dataset is available.
417 catalogs : list of lsst.afw.table.SourceCatalog
418 Catalogs to combine into per-visit catalogs.
419 visits : list of visit id (usually int)
420 list of visit identifiers, one-to-one correspondent with catalogs.
425 dict of visit: catalog of all sources from all CCDs of that visit.
428 for v, cat
in zip(visits, catalogs):
429 visit_dict[v].extend(cat)
432 visit_dict[v] = visit_dict[v].copy(deep=
True)
438 faint, bright, old_PA1, new_PA1,
439 name='', outdir='.plots'):
440 """Plot various distributions of fluxes and magnitudes.
444 plt : matplotlib.pyplot instance
445 pyplot instance to plot with
450 old_weighted_rms : np.array
451 old rms weighted by the mean (rms(data)/mean(data))
452 new_weighted_rms : np.array
453 old rms weighted by the mean (rms(data)/mean(data))
455 Faint end of range that PA1 was computed from.
457 Bright end of range that PA1 was computed from.
459 Old value of PA1, to plot as horizontal line.
461 New value of PA1, to plot as horizontal line.
463 Name to include in plot titles and save files.
464 outdir : str, optional
465 Directory to write the saved plots to.
469 seaborn.set_style(
'whitegrid')
475 plt.plot(old_mag, old_weighted_rms,
'.', color=old_color, label=
'old')
476 plt.plot(new_mag, new_weighted_rms,
'.', color=new_color, label=
'new')
477 plt.axvline(faint, ls=
':', color=old_color)
478 plt.axvline(bright, ls=
':', color=old_color)
479 plt.axhline(old_PA1, ls=
'--', color=old_color)
480 plt.axhline(new_PA1, ls=
'--', color=new_color)
481 plt.legend(loc=
'upper left')
482 plt.title(
'Where is the systematic flux rms limit?')
483 plt.xlabel(
'magnitude')
484 plt.ylabel(
'rms/mean per source')
485 filename = os.path.join(outdir,
'{}-photometry-PA1.pdf')
486 plt.savefig(filename.format(name))
489 seaborn.distplot(old_weighted_rms, fit=scipy.stats.lognorm, kde=
False, label=
"old", color=old_color)
490 seaborn.distplot(new_weighted_rms, fit=scipy.stats.lognorm, kde=
False, label=
"new", color=new_color)
491 plt.title(
'Source RMS pre/post-jointcal')
492 plt.xlabel(
'rms(flux)/mean(flux)')
494 plt.legend(loc=
'upper right')
495 filename = os.path.join(outdir,
'{}-photometry-rms.pdf')
496 plt.savefig(filename.format(name))
500 name='', outdir='.plots'):
502 Various plots of the difference between old and new Wcs.
506 plt : matplotlib.pyplot instance
507 pyplot instance to plot with.
508 data_refs : list of lsst.daf.persistence.butlerSubset.ButlerDataRef
509 A list of data refs to plot.
510 visits : list of visit id (usually int)
511 list of visit identifiers, one-to-one correspondent with catalogs.
512 old_wcs_list : list of lsst.afw.image.wcs.Wcs
513 A list of the old (pre-jointcal) WCSs, one-to-one corresponding to data_refs.
514 per_ccd_plot : bool, optional
515 Make per-ccd plots of the "wcs different" (warning: slow!)
517 Name to include in plot titles and save files.
518 outdir : str, optional
519 Directory to write the saved plots to.
526 for i, ref
in enumerate(data_refs):
527 md = ref.get(
'calexp_md')
529 plot_wcs(plt, old_wcs_list[i], ref.get(
'jointcal_wcs'),
530 dims.getX(), dims.getY(),
531 center=(md.getScalar(
'CRVAL1'), md.getScalar(
'CRVAL2')), name=
'dataRef %d'%i,
536 """Return num x/y grid coordinates for wcs1 and wcs2."""
537 x = np.linspace(0, x_dim, num)
538 y = np.linspace(0, y_dim, num)
541 return x1, y1, x2, y2
545 """Convert two arrays of x/y points into an on-sky grid."""
546 xout = np.zeros((xv.shape[0], yv.shape[0]))
547 yout = np.zeros((xv.shape[0], yv.shape[0]))
548 for i, x
in enumerate(xv):
549 for j, y
in enumerate(yv):
550 sky = wcs.pixelToSky(x, y)
551 xout[i, j] = sky.getRa()
552 yout[i, j] = sky.getDec()
558 Make quiver plots of the WCS deltas for each CCD in each visit.
562 plt : matplotlib.pyplot instance
563 pyplot instance to plot with.
564 data_refs : list of lsst.daf.persistence.butlerSubset.ButlerDataRef
565 A list of data refs to plot.
566 visits : list of visit id (usually int)
567 list of visit identifiers, one-to-one correspondent with catalogs.
568 old_wcs_list : list of lsst.afw.image.wcs.Wcs
569 A list of the old (pre-jointcal) WCSs, one-to-one corresponding to data_refs.
571 Name to include in plot titles and save files.
572 outdir : str, optional
573 Directory to write the saved plots to.
579 ax = fig.add_subplot(111)
580 for old_wcs, ref
in zip(old_wcs_list, data_refs):
581 if ref.dataId[
'visit'] != visit:
583 md = ref.get(
'calexp_md')
586 dims.getX(), dims.getY())
589 length = (0.1*u.arcsecond).
to(u.radian).value
590 ax.quiverkey(Q, 0.9, 0.95, length,
'0.1 arcsec', coordinates=
'figure', labelpos=
'W')
593 plt.title(
'visit: {}'.
format(visit))
594 filename = os.path.join(outdir,
'{}-{}-quivers.pdf')
595 plt.savefig(filename.format(name, visit))
600 Plot the delta between wcs1 and wcs2 as vector arrows.
605 Matplotlib axis instance to plot to.
606 wcs1 : lsst.afw.image.wcs.Wcs
607 First WCS to compare.
608 wcs2 : lsst.afw.image.wcs.Wcs
609 Second WCS to compare.
611 Size of array in X-coordinate to make the grid over.
613 Size of array in Y-coordinate to make the grid over.
619 return ax.quiver(x1, y1, uu, vv, units=
'x', pivot=
'tail', scale=1e-3, width=1e-5)
623 """Plot the magnitude of the WCS change between old and new visits as a heat map.
627 plt : matplotlib.pyplot instance
628 pyplot instance to plot with.
629 data_refs : list of lsst.daf.persistence.butlerSubset.ButlerDataRef
630 A list of data refs to plot.
631 visits : list of visit id (usually int)
632 list of visit identifiers, one-to-one correspondent with catalogs.
633 old_wcs_list : list of lsst.afw.image.wcs.Wcs
634 A list of the old (pre-jointcal) WCSs, one-to-one corresponding to data_refs.
636 Name to include in plot titles and save files.
637 outdir : str, optional
638 Directory to write the saved plots to.
642 fig.set_tight_layout(
True)
643 ax = fig.add_subplot(111)
649 for old_wcs, ref
in zip(old_wcs_list, data_refs):
650 if ref.dataId[
'visit'] != visit:
652 md = ref.get(
'calexp_md')
655 old_wcs, ref.get(
'jointcal_wcs'))
658 extent = (x1[0, 0], x1[-1, -1], y1[0, 0], y1[-1, -1])
659 xmin =
min(x1.min(), xmin)
660 ymin =
min(y1.min(), ymin)
661 xmax =
max(x1.max(), xmax)
662 ymax =
max(y1.max(), ymax)
663 magnitude = (np.linalg.norm((uu, vv), axis=0)*u.radian).
to(u.arcsecond).value
664 img = ax.imshow(magnitude, vmin=0, vmax=0.3,
665 aspect=
'auto', extent=extent, cmap=plt.get_cmap(
'magma'))
671 cbar = plt.colorbar(img)
672 cbar.ax.set_ylabel(
'distortion (arcseconds)')
677 plt.title(
'visit: {}'.
format(visit))
678 filename = os.path.join(outdir,
'{}-{}-heatmap.pdf')
679 plt.savefig(filename.format(name, visit))
682 def plot_wcs(plt, wcs1, wcs2, x_dim, y_dim, center=(0, 0), name=
"", outdir=
'.plots'):
683 """Plot the "distortion map": wcs1-wcs2 delta of points in the CCD grid.
687 plt : matplotlib.pyplot instance
688 pyplot instance to plot with.
689 wcs1 : lsst.afw.image.wcs.Wcs
690 First WCS to compare.
691 wcs2 : lsst.afw.image.wcs.Wcs
692 Second WCS to compare.
694 Size of array in X-coordinate to make the grid over.
696 Size of array in Y-coordinate to make the grid over.
697 center : tuple, optional
698 Center of the data, in on-chip coordinates.
700 Name to include in plot titles and save files.
701 outdir : str, optional
702 Directory to write the saved plots to.
708 plt.plot((x1 - x2) + center[0], (y1 - y2) + center[1],
'-')
709 plt.xlabel(
'delta RA (arcsec)')
710 plt.ylabel(
'delta Dec (arcsec)')
712 filename = os.path.join(outdir,
'{}-wcs.pdf')
713 plt.savefig(filename.format(name))
717 new_rms_relative, new_rms_absolute,
718 old_rel_total, old_abs_total, new_rel_total, new_abs_total,
719 name="", outdir='.plots'):
720 """Plot histograms of the source separations and their RMS values.
724 plt : matplotlib.pyplot instance
725 pyplot instance to plot with.
726 old_rms_relative : np.array
727 old relative rms/star
728 old_rms_absolute : np.array
729 old absolute rms/star
730 new_rms_relative : np.array
731 new relative rms/star
732 new_rms_absolute : np.array
733 new absolute rms/star
734 old_rel_total : float
735 old relative rms over all stars
736 old_abs_total : float
737 old absolute rms over all stars
738 new_rel_total : float
739 new relative rms over all stars
740 new_abs_total : float
741 new absolute rms over all stars
743 Name to include in plot titles and save files.
744 outdir : str, optional
745 Directory to write the saved plots to.
753 plotOptions = {
'lw': 2,
'range': (0, 0.1)*u.arcsecond,
'normed':
True,
754 'bins': 30,
'histtype':
'step'}
756 plt.title(
'relative vs. absolute: %d vs. %d'%(len(old_rms_relative), len(old_rms_absolute)))
758 plt.hist(old_rms_absolute, color=color_abs, ls=ls_old, label=
'old abs', **plotOptions)
759 plt.hist(new_rms_absolute, color=color_abs, ls=ls_new, label=
'new abs', **plotOptions)
761 plt.hist(old_rms_relative, color=color_rel, ls=ls_old, label=
'old rel', **plotOptions)
762 plt.hist(new_rms_relative, color=color_rel, ls=ls_new, label=
'new rel', **plotOptions)
764 plt.axvline(x=old_abs_total.value, linewidth=1.5, color=color_abs, ls=ls_old)
765 plt.axvline(x=new_abs_total.value, linewidth=1.5, color=color_abs, ls=ls_new)
766 plt.axvline(x=old_rel_total.value, linewidth=1.5, color=color_rel, ls=ls_old)
767 plt.axvline(x=new_rel_total.value, linewidth=1.5, color=color_rel, ls=ls_new)
769 plt.xlim(plotOptions[
'range'])
770 plt.xlabel(
'arcseconds')
771 plt.legend(loc=
'best')
772 filename = os.path.join(outdir,
'{}-histogram.pdf')
773 plt.savefig(filename.format(name))
def _make_match_dict(self, reference, visit_catalogs, photoCalibs, refcalib=None)
def compute_rms(self, data_refs, reference)
def _photometric_rms(self, sn_cut=300, magnitude_range=3)
def make_plots(self, data_refs, old_wcs_list, name='', interactive=False, per_ccd_plot=False, outdir='.plots')
def _make_visit_catalogs(self, catalogs, visits)
def __init__(self, match_radius=0.1 *arcseconds, flux_limit=100.0, do_photometry=True, do_astrometry=True, verbose=False)
static Log getLogger(Log const &logger)
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
lsst::geom::Box2I bboxFromMetadata(daf::base::PropertySet &metadata)
Determine the image bounding box from its metadata (FITS header)
void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList)
Update sky coordinates in a collection of source objects.
std::vector< Match< typename Cat1::Record, typename Cat2::Record > > matchRaDec(Cat1 const &cat1, Cat2 const &cat2, lsst::geom::Angle radius, MatchControl const &mc=MatchControl())
Compute all tuples (s1,s2,d) where s1 belings to cat1, s2 belongs to cat2 and d, the distance between...
def plot_all_wcs_deltas(plt, data_refs, visits, old_wcs_list, per_ccd_plot=False, name='', outdir='.plots')
def plot_wcs_quivers(ax, wcs1, wcs2, x_dim, y_dim)
def plot_wcs(plt, wcs1, wcs2, x_dim, y_dim, center=(0, 0), name="", outdir='.plots')
def wcs_convert(xv, yv, wcs)
def make_xy_wcs_grid(x_dim, y_dim, wcs1, wcs2, num=50)
def plot_flux_distributions(plt, old_mag, new_mag, old_weighted_rms, new_weighted_rms, faint, bright, old_PA1, new_PA1, name='', outdir='.plots')
def plot_rms_histogram(plt, old_rms_relative, old_rms_absolute, new_rms_relative, new_rms_absolute, old_rel_total, old_abs_total, new_rel_total, new_abs_total, name="", outdir='.plots')
def plot_all_wcs_quivers(plt, data_refs, visits, old_wcs_list, name, outdir='.plots')
def plot_wcs_magnitude(plt, data_refs, visits, old_wcs_list, name, outdir='.plots')
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)