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'])
51 Compute statistics on jointcal-processed data, and optionally generate plots.
55 Instantiate JointcalStatistics and call compute_rms() to get the relevant
56 statistics for e.g. unittests, and call make_plots() to generate a suite of
60 def __init__(self, match_radius=0.1*arcseconds, flux_limit=100.0,
61 do_photometry=True, do_astrometry=True,
66 match_radius : lsst.geom.Angle
67 match sources within this radius for RMS statistics
69 Signal/Noise (flux/fluxErr) for sources to be included in the RMS cross-match.
70 100 is a balance between good centroids and enough sources.
71 do_photometry : bool, optional
72 Perform calculations/make plots for photometric metrics.
73 do_astrometry : bool, optional
74 Perform calculations/make plots for astrometric metrics.
75 verbose : bool, optional
87 Match all data_refs to compute the RMS, for all detections above self.flux_limit.
91 data_refs : list of lsst.daf.persistence.butlerSubset.ButlerDataRef
92 A list of data refs to do the calculations between.
93 reference : lsst reference catalog
94 reference catalog to do absolute matching against.
100 Post-jointcal relative RMS of the matched sources.
102 Post-jointcal absolute RMS of matched sources.
104 Post-jointcal photometric repeatability (PA1 from the SRD).
108 self.
filters = [ref.get(
'calexp_filter').getName()
for ref
in data_refs]
111 def compute(catalogs, photoCalibs):
112 """Compute the relative and absolute matches in distance and flux."""
122 refcalib = photoCalibs[0]
if photoCalibs != []
else None
123 dist_rel, flux_rel, ref_flux_rel, source_rel = self.
_make_match_dict(refcat,
127 dist_abs, flux_abs, ref_flux_abs, source_abs = self.
_make_match_dict(reference,
132 ref_flux =
MatchDict(ref_flux_rel, ref_flux_abs)
133 source =
MatchDict(source_rel, source_abs)
134 return dist, flux, ref_flux, source
136 old_cats = [ref.get(
'src')
for ref
in data_refs]
141 old_calibs = [ref.get(
'calexp_photoCalib')
for ref
in data_refs]
143 self.old_dist, self.old_flux, self.old_ref_flux, self.
old_source = compute(old_cats, old_calibs)
146 new_cats = [ref.get(
'src')
for ref
in data_refs]
149 new_wcss = [ref.get(
'jointcal_wcs')
for ref
in data_refs]
152 new_calibs = [ref.get(
'jointcal_photoCalib')
for ref
in data_refs]
154 for wcs, cat
in zip(new_wcss, new_cats):
158 self.new_dist, self.new_flux, self.new_ref_flux, self.
new_source = compute(new_cats, new_calibs)
161 print(
'old, new relative distance matches:',
162 len(self.old_dist.relative), len(self.new_dist.relative))
163 print(
'old, new absolute distance matches:',
164 len(self.old_dist.absolute), len(self.new_dist.absolute))
165 print(
'old, new relative flux matches:',
166 len(self.old_flux.relative), len(self.new_flux.relative))
167 print(
'old, new absolute flux matches:',
168 len(self.old_flux.absolute), len(self.new_flux.absolute))
176 """Compute the total rms across all sources."""
177 total = sum(sum(dd**2)
for dd
in data.values())
178 n = sum(len(dd)
for dd
in data.values())
179 return np.sqrt(total/n)
190 print(
"relative, absolute astrometry:",
194 print(
"Photometric Accuracy (PA1):", self.
new_PA1)
195 Rms_result = collections.namedtuple(
"Rms_result", [
"dist_relative",
"dist_absolute",
"pa1"])
199 name='', interactive=False, per_ccd_plot=False, outdir='.plots'):
201 Make plots of various quantites to help with debugging.
202 Requires that `compute_rms()` was run first.
206 data_refs : list of lsst.daf.persistence.butlerSubset.ButlerDataRef
207 A list of data refs to do the calculations between.
208 old_wcs_list : list of lsst.afw.image.wcs.Wcs
209 A list of the old (pre-jointcal) WCSs, one-to-one corresponding to data_refs.
211 Name to include in plot titles and save files.
213 Turn on matplotlib interactive mode and drop into a debugger when
214 plotting is finished. Otherwise, use a non-interactive backend.
216 Plot the WCS per CCD (takes longer and generates many plots for a large camera)
218 directory to save plots to.
224 matplotlib.use(
'pdf')
226 import matplotlib.pyplot
as plt
227 import astropy.visualization
229 astropy.visualization.quantity_support()
233 self.
log.
info(
"N data_refs: %d", len(data_refs))
239 name=name, outdir=outdir)
243 def rms_per_source(data):
244 """Each element of data must already be the "delta" of whatever measurement."""
245 return (np.sqrt([np.mean(dd**2)
for dd
in data.values()])*u.radian).
to(u.arcsecond)
256 new_dist_rms.relative, new_dist_rms.absolute,
259 name=name, outdir=outdir)
262 per_ccd_plot=per_ccd_plot,
263 name=name, outdir=outdir)
270 def _photometric_rms(self, sn_cut=300, magnitude_range=3):
272 Compute the photometric RMS and the photometric repeatablity values (PA1).
277 The minimum signal/noise for sources to be included in the PA1 calculation.
278 magnitude_range : float
279 The range of magnitudes above sn_cut to include in the PA1 calculation.
281 def rms(flux, ref_flux):
282 return np.sqrt([np.mean((ref_flux[dd] - flux[dd])**2)
for dd
in flux])
288 self.
old_ref = np.fromiter(self.old_ref_flux.absolute.values(), dtype=float)
289 self.
new_ref = np.fromiter(self.new_ref_flux.absolute.values(), dtype=float)
293 def signal_to_noise(sources, flux_key='slot_PsfFlux_instFlux', sigma_key='slot_PsfFlux_instFluxErr'):
294 """Compute the mean signal/noise per source from a MatchDict of SourceRecords."""
295 result = np.empty(len(sources))
296 for i, src
in enumerate(sources.values()):
297 result[i] = np.mean([x[flux_key]/x[sigma_key]
for x
in src])
300 old_sn = signal_to_noise(self.
old_source.absolute)
313 def _make_match_dict(self, reference, visit_catalogs, photoCalibs, refcalib=None):
315 Return several dicts of sourceID:[values] over the catalogs, to be used in RMS calculations.
319 reference : lsst.afw.table.SourceCatalog
320 Catalog to do the matching against.
321 visit_catalogs : list of lsst.afw.table.SourceCatalog
322 Visit source catalogs (values() produced by _make_visit_catalogs)
323 to cross-match against reference.
324 photoCalibs : list of lsst.afw.image.PhotoCalib
325 Exposure PhotoCalibs, 1-1 coorespondent with visit_catalogs.
326 refcalib : lsst.afw.image.PhotoCalib or None
327 Pass a PhotoCalib here to use it to compute nanojansky from the
328 reference catalog ADU slot_flux.
333 dict of sourceID: array(separation distances for that source)
335 dict of sourceID: array(fluxes (nJy) for that source)
337 dict of sourceID: flux (nJy) of the reference object
339 dict of sourceID: list(each SourceRecord that was position-matched
343 if photoCalibs == []:
344 photoCalibs = [[]]*len(visit_catalogs)
346 distances = collections.defaultdict(list)
347 fluxes = collections.defaultdict(list)
349 sources = collections.defaultdict(list)
350 if 'slot_CalibFlux_instFlux' in reference.schema:
351 ref_flux_key =
'slot_CalibFlux'
353 ref_flux_key =
'{}_flux'
355 def get_fluxes(photoCalib, match):
356 """Return (flux, ref_flux) or None if either is invalid."""
358 flux = match[1][
'slot_CalibFlux_instFlux']
362 flux = photoCalib.instFluxToNanojansky(match[1],
"slot_CalibFlux").value
365 if 'slot' in ref_flux_key:
366 ref_flux = match[0][ref_flux_key+
'_instFlux']
370 ref_flux = refcalib.instFluxToNanojansky(match[0], ref_flux_key).value
373 ref_flux = match[0][ref_flux_key.format(filt)]
377 Flux = collections.namedtuple(
'Flux', (
'flux',
'ref_flux'))
378 return Flux(flux, ref_flux)
380 for cat, photoCalib, filt
in zip(visit_catalogs, photoCalibs, self.
filters):
381 good = (cat.get(
'base_PsfFlux_instFlux')/cat.get(
'base_PsfFlux_instFluxErr')) > self.
flux_limit
383 good &= (cat.get(
'base_ClassificationExtendedness_value') == 0)
387 flux = get_fluxes(photoCalib, m)
391 fluxes[m[0].getId()].
append(flux.flux)
393 ref_fluxes[m[0].getId()] = flux.ref_flux
397 distances[m[0].getId()].
append(m[2])
399 sources[m[0].getId()].
append(m[1])
401 for source
in distances:
402 distances[source] = np.array(distances[source])
403 for source
in fluxes:
404 fluxes[source] = np.array(fluxes[source])
406 return distances, fluxes, ref_fluxes, sources
408 def _make_visit_catalogs(self, catalogs, visits):
410 Merge all catalogs from the each visit.
411 NOTE: creating this structure is somewhat slow, and will be unnecessary
412 once a full-visit composite dataset is available.
416 catalogs : list of lsst.afw.table.SourceCatalog
417 Catalogs to combine into per-visit catalogs.
418 visits : list of visit id (usually int)
419 list of visit identifiers, one-to-one correspondent with catalogs.
424 dict of visit: catalog of all sources from all CCDs of that visit.
427 for v, cat
in zip(visits, catalogs):
428 visit_dict[v].extend(cat)
431 visit_dict[v] = visit_dict[v].copy(deep=
True)
437 faint, bright, old_PA1, new_PA1,
438 name='', outdir='.plots'):
439 """Plot various distributions of fluxes and magnitudes.
443 plt : matplotlib.pyplot instance
444 pyplot instance to plot with
449 old_weighted_rms : np.array
450 old rms weighted by the mean (rms(data)/mean(data))
451 new_weighted_rms : np.array
452 old rms weighted by the mean (rms(data)/mean(data))
454 Faint end of range that PA1 was computed from.
456 Bright end of range that PA1 was computed from.
458 Old value of PA1, to plot as horizontal line.
460 New value of PA1, to plot as horizontal line.
462 Name to include in plot titles and save files.
463 outdir : str, optional
464 Directory to write the saved plots to.
468 seaborn.set_style(
'whitegrid')
474 plt.plot(old_mag, old_weighted_rms,
'.', color=old_color, label=
'old')
475 plt.plot(new_mag, new_weighted_rms,
'.', color=new_color, label=
'new')
476 plt.axvline(faint, ls=
':', color=old_color)
477 plt.axvline(bright, ls=
':', color=old_color)
478 plt.axhline(old_PA1, ls=
'--', color=old_color)
479 plt.axhline(new_PA1, ls=
'--', color=new_color)
480 plt.legend(loc=
'upper left')
481 plt.title(
'Where is the systematic flux rms limit?')
482 plt.xlabel(
'magnitude')
483 plt.ylabel(
'rms/mean per source')
484 filename = os.path.join(outdir,
'{}-photometry-PA1.pdf')
485 plt.savefig(filename.format(name))
488 seaborn.distplot(old_weighted_rms, fit=scipy.stats.lognorm, kde=
False, label=
"old", color=old_color)
489 seaborn.distplot(new_weighted_rms, fit=scipy.stats.lognorm, kde=
False, label=
"new", color=new_color)
490 plt.title(
'Source RMS pre/post-jointcal')
491 plt.xlabel(
'rms(flux)/mean(flux)')
493 plt.legend(loc=
'upper right')
494 filename = os.path.join(outdir,
'{}-photometry-rms.pdf')
495 plt.savefig(filename.format(name))
499 name='', outdir='.plots'):
501 Various plots of the difference between old and new Wcs.
505 plt : matplotlib.pyplot instance
506 pyplot instance to plot with.
507 data_refs : list of lsst.daf.persistence.butlerSubset.ButlerDataRef
508 A list of data refs to plot.
509 visits : list of visit id (usually int)
510 list of visit identifiers, one-to-one correspondent with catalogs.
511 old_wcs_list : list of lsst.afw.image.wcs.Wcs
512 A list of the old (pre-jointcal) WCSs, one-to-one corresponding to data_refs.
513 per_ccd_plot : bool, optional
514 Make per-ccd plots of the "wcs different" (warning: slow!)
516 Name to include in plot titles and save files.
517 outdir : str, optional
518 Directory to write the saved plots to.
525 for i, ref
in enumerate(data_refs):
526 md = ref.get(
'calexp_md')
528 plot_wcs(plt, old_wcs_list[i], ref.get(
'jointcal_wcs'),
529 dims.getX(), dims.getY(),
530 center=(md.getScalar(
'CRVAL1'), md.getScalar(
'CRVAL2')), name=
'dataRef %d'%i,
535 """Return num x/y grid coordinates for wcs1 and wcs2."""
536 x = np.linspace(0, x_dim, num)
537 y = np.linspace(0, y_dim, num)
540 return x1, y1, x2, y2
544 """Convert two arrays of x/y points into an on-sky grid."""
545 xout = np.zeros((xv.shape[0], yv.shape[0]))
546 yout = np.zeros((xv.shape[0], yv.shape[0]))
547 for i, x
in enumerate(xv):
548 for j, y
in enumerate(yv):
549 sky = wcs.pixelToSky(x, y)
550 xout[i, j] = sky.getRa()
551 yout[i, j] = sky.getDec()
557 Make quiver plots of the WCS deltas for each CCD in each visit.
561 plt : matplotlib.pyplot instance
562 pyplot instance to plot with.
563 data_refs : list of lsst.daf.persistence.butlerSubset.ButlerDataRef
564 A list of data refs to plot.
565 visits : list of visit id (usually int)
566 list of visit identifiers, one-to-one correspondent with catalogs.
567 old_wcs_list : list of lsst.afw.image.wcs.Wcs
568 A list of the old (pre-jointcal) WCSs, one-to-one corresponding to data_refs.
570 Name to include in plot titles and save files.
571 outdir : str, optional
572 Directory to write the saved plots to.
578 ax = fig.add_subplot(111)
579 for old_wcs, ref
in zip(old_wcs_list, data_refs):
580 if ref.dataId[
'visit'] != visit:
582 md = ref.get(
'calexp_md')
585 dims.getX(), dims.getY())
588 length = (0.1*u.arcsecond).
to(u.radian).value
589 ax.quiverkey(Q, 0.9, 0.95, length,
'0.1 arcsec', coordinates=
'figure', labelpos=
'W')
592 plt.title(
'visit: {}'.
format(visit))
593 filename = os.path.join(outdir,
'{}-{}-quivers.pdf')
594 plt.savefig(filename.format(name, visit))
599 Plot the delta between wcs1 and wcs2 as vector arrows.
604 Matplotlib axis instance to plot to.
605 wcs1 : lsst.afw.image.wcs.Wcs
606 First WCS to compare.
607 wcs2 : lsst.afw.image.wcs.Wcs
608 Second WCS to compare.
610 Size of array in X-coordinate to make the grid over.
612 Size of array in Y-coordinate to make the grid over.
618 return ax.quiver(x1, y1, uu, vv, units=
'x', pivot=
'tail', scale=1e-3, width=1e-5)
622 """Plot the magnitude of the WCS change between old and new visits as a heat map.
626 plt : matplotlib.pyplot instance
627 pyplot instance to plot with.
628 data_refs : list of lsst.daf.persistence.butlerSubset.ButlerDataRef
629 A list of data refs to plot.
630 visits : list of visit id (usually int)
631 list of visit identifiers, one-to-one correspondent with catalogs.
632 old_wcs_list : list of lsst.afw.image.wcs.Wcs
633 A list of the old (pre-jointcal) WCSs, one-to-one corresponding to data_refs.
635 Name to include in plot titles and save files.
636 outdir : str, optional
637 Directory to write the saved plots to.
641 fig.set_tight_layout(
True)
642 ax = fig.add_subplot(111)
648 for old_wcs, ref
in zip(old_wcs_list, data_refs):
649 if ref.dataId[
'visit'] != visit:
651 md = ref.get(
'calexp_md')
654 old_wcs, ref.get(
'jointcal_wcs'))
657 extent = (x1[0, 0], x1[-1, -1], y1[0, 0], y1[-1, -1])
658 xmin =
min(x1.min(), xmin)
659 ymin =
min(y1.min(), ymin)
660 xmax =
max(x1.max(), xmax)
661 ymax =
max(y1.max(), ymax)
662 magnitude = (np.linalg.norm((uu, vv), axis=0)*u.radian).
to(u.arcsecond).value
663 img = ax.imshow(magnitude, vmin=0, vmax=0.3,
664 aspect=
'auto', extent=extent, cmap=plt.get_cmap(
'magma'))
670 cbar = plt.colorbar(img)
671 cbar.ax.set_ylabel(
'distortion (arcseconds)')
676 plt.title(
'visit: {}'.
format(visit))
677 filename = os.path.join(outdir,
'{}-{}-heatmap.pdf')
678 plt.savefig(filename.format(name, visit))
681 def plot_wcs(plt, wcs1, wcs2, x_dim, y_dim, center=(0, 0), name=
"", outdir=
'.plots'):
682 """Plot the "distortion map": wcs1-wcs2 delta of points in the CCD grid.
686 plt : matplotlib.pyplot instance
687 pyplot instance to plot with.
688 wcs1 : lsst.afw.image.wcs.Wcs
689 First WCS to compare.
690 wcs2 : lsst.afw.image.wcs.Wcs
691 Second WCS to compare.
693 Size of array in X-coordinate to make the grid over.
695 Size of array in Y-coordinate to make the grid over.
696 center : tuple, optional
697 Center of the data, in on-chip coordinates.
699 Name to include in plot titles and save files.
700 outdir : str, optional
701 Directory to write the saved plots to.
707 plt.plot((x1 - x2) + center[0], (y1 - y2) + center[1],
'-')
708 plt.xlabel(
'delta RA (arcsec)')
709 plt.ylabel(
'delta Dec (arcsec)')
711 filename = os.path.join(outdir,
'{}-wcs.pdf')
712 plt.savefig(filename.format(name))
716 new_rms_relative, new_rms_absolute,
717 old_rel_total, old_abs_total, new_rel_total, new_abs_total,
718 name="", outdir='.plots'):
719 """Plot histograms of the source separations and their RMS values.
723 plt : matplotlib.pyplot instance
724 pyplot instance to plot with.
725 old_rms_relative : np.array
726 old relative rms/star
727 old_rms_absolute : np.array
728 old absolute rms/star
729 new_rms_relative : np.array
730 new relative rms/star
731 new_rms_absolute : np.array
732 new absolute rms/star
733 old_rel_total : float
734 old relative rms over all stars
735 old_abs_total : float
736 old absolute rms over all stars
737 new_rel_total : float
738 new relative rms over all stars
739 new_abs_total : float
740 new absolute rms over all stars
742 Name to include in plot titles and save files.
743 outdir : str, optional
744 Directory to write the saved plots to.
752 plotOptions = {
'lw': 2,
'range': (0, 0.1)*u.arcsecond,
'normed':
True,
753 'bins': 30,
'histtype':
'step'}
755 plt.title(
'relative vs. absolute: %d vs. %d'%(len(old_rms_relative), len(old_rms_absolute)))
757 plt.hist(old_rms_absolute, color=color_abs, ls=ls_old, label=
'old abs', **plotOptions)
758 plt.hist(new_rms_absolute, color=color_abs, ls=ls_new, label=
'new abs', **plotOptions)
760 plt.hist(old_rms_relative, color=color_rel, ls=ls_old, label=
'old rel', **plotOptions)
761 plt.hist(new_rms_relative, color=color_rel, ls=ls_new, label=
'new rel', **plotOptions)
763 plt.axvline(x=old_abs_total.value, linewidth=1.5, color=color_abs, ls=ls_old)
764 plt.axvline(x=new_abs_total.value, linewidth=1.5, color=color_abs, ls=ls_new)
765 plt.axvline(x=old_rel_total.value, linewidth=1.5, color=color_rel, ls=ls_old)
766 plt.axvline(x=new_rel_total.value, linewidth=1.5, color=color_rel, ls=ls_new)
768 plt.xlim(plotOptions[
'range'])
769 plt.xlabel(
'arcseconds')
770 plt.legend(loc=
'best')
771 filename = os.path.join(outdir,
'{}-histogram.pdf')
772 plt.savefig(filename.format(name))