LSSTApplications  19.0.0-10-g920eed2,19.0.0-11-g48a0200+2,19.0.0-18-gfc4e62b+13,19.0.0-2-g3b2f90d+2,19.0.0-2-gd671419+5,19.0.0-20-g5a5a17ab+11,19.0.0-21-g2644856+13,19.0.0-23-g84eeccb+1,19.0.0-24-g878c510+1,19.0.0-25-g6c8df7140,19.0.0-25-gb330496+1,19.0.0-3-g2b32d65+5,19.0.0-3-g8227491+12,19.0.0-3-g9c54d0d+12,19.0.0-3-gca68e65+8,19.0.0-3-gcfc5f51+5,19.0.0-3-ge110943+11,19.0.0-3-ge74d124,19.0.0-3-gfe04aa6+13,19.0.0-30-g9c3fd16+1,19.0.0-4-g06f5963+5,19.0.0-4-g3d16501+13,19.0.0-4-g4a9c019+5,19.0.0-4-g5a8b323,19.0.0-4-g66397f0+1,19.0.0-4-g8278b9b+1,19.0.0-4-g8557e14,19.0.0-4-g8964aba+13,19.0.0-4-ge404a01+12,19.0.0-5-g40f3a5a,19.0.0-5-g4db63b3,19.0.0-5-gfb03ce7+13,19.0.0-6-gbaebbfb+12,19.0.0-61-gec4c6e08+1,19.0.0-7-g039c0b5+11,19.0.0-7-gbea9075+4,19.0.0-7-gc567de5+13,19.0.0-71-g41c0270,19.0.0-9-g2f02add+1,19.0.0-9-g463f923+12,w.2020.22
LSSTDataManagementBasePackage
utils.py
Go to the documentation of this file.
1 # This file is part of jointcal.
2 #
3 # Developed for the LSST Data Management System.
4 # This product includes software developed by the LSST Project
5 # (https://www.lsst.org).
6 # See the COPYRIGHT file at the top-level directory of this distribution
7 # for details of code ownership.
8 #
9 # This program is free software: you can redistribute it and/or modify
10 # it under the terms of the GNU General Public License as published by
11 # the Free Software Foundation, either version 3 of the License, or
12 # (at your option) any later version.
13 #
14 # This program is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU General Public License for more details.
18 #
19 # You should have received a copy of the GNU General Public License
20 # along with this program. If not, see <https://www.gnu.org/licenses/>.
21 
22 """
23 Statistics of jointcal vs. single-frame procesing and diagnostic plots.
24 
25 Notes
26 -----
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:
29 
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
32 """
33 import collections
34 import os
35 
36 import numpy as np
37 from astropy import units as u
38 
39 import lsst.log
40 import lsst.afw.table
41 import lsst.afw.image
42 from lsst.geom import arcseconds
43 
44 __all__ = ['JointcalStatistics']
45 
46 MatchDict = collections.namedtuple('MatchDict', ['relative', 'absolute'])
47 
48 
50  """
51  Compute statistics on jointcal-processed data, and optionally generate plots.
52 
53  Notes
54  -----
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
57  diagnostic plots.
58  """
59 
60  def __init__(self, match_radius=0.1*arcseconds, flux_limit=100.0,
61  do_photometry=True, do_astrometry=True,
62  verbose=False):
63  """
64  Parameters
65  ----------
66  match_radius : lsst.geom.Angle
67  match sources within this radius for RMS statistics
68  flux_limit : float
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
76  Print extra things
77  """
78  self.match_radius = match_radius
79  self.flux_limit = flux_limit
80  self.do_photometry = do_photometry
81  self.do_astrometry = do_astrometry
82  self.verbose = verbose
83  self.log = lsst.log.Log.getLogger('JointcalStatistics')
84 
85  def compute_rms(self, data_refs, reference):
86  """
87  Match all data_refs to compute the RMS, for all detections above self.flux_limit.
88 
89  Parameters
90  ----------
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.
95 
96  Returns
97  -------
98  namedtuple:
99  astropy.Quantity
100  Post-jointcal relative RMS of the matched sources.
101  astropy.Quantity
102  Post-jointcal absolute RMS of matched sources.
103  float
104  Post-jointcal photometric repeatability (PA1 from the SRD).
105  """
106 
107  # DECAM doesn't have "filter" in its registry, so we have to get the filter names directly.
108  self.filters = [ref.get('calexp_filter').getName() for ref in data_refs]
109  self.visits_per_dataRef = [ref.dataId['visit'] for ref in data_refs]
110 
111  def compute(catalogs, photoCalibs):
112  """Compute the relative and absolute matches in distance and flux."""
113  visit_catalogs = self._make_visit_catalogs(catalogs, self.visits_per_dataRef)
114  catalogs = [visit_catalogs[x] for x in self.visits_per_dataRef]
115  # use the first catalog as the relative reference catalog
116  # NOTE: The "first" catalog depends on the original ordering of the data_refs.
117  # NOTE: Thus, because I'm doing a many-1 match in _make_match_dict,
118  # the number of matches (and thus the details of the match statistics)
119  # will change if the data_refs are ordered differently.
120  # All the more reason to use a proper n-way matcher here. See: DM-8664
121  refcat = catalogs[0]
122  refcalib = photoCalibs[0] if photoCalibs != [] else None
123  dist_rel, flux_rel, ref_flux_rel, source_rel = self._make_match_dict(refcat,
124  catalogs[1:],
125  photoCalibs[1:],
126  refcalib=refcalib)
127  dist_abs, flux_abs, ref_flux_abs, source_abs = self._make_match_dict(reference,
128  catalogs,
129  photoCalibs)
130  dist = MatchDict(dist_rel, dist_abs)
131  flux = MatchDict(flux_rel, flux_abs)
132  ref_flux = MatchDict(ref_flux_rel, ref_flux_abs)
133  source = MatchDict(source_rel, source_abs)
134  return dist, flux, ref_flux, source
135 
136  old_cats = [ref.get('src') for ref in data_refs]
137  # NOTE: build photoCalibs from existing old Calib objects.
138  # TODO: we can make this a listcomp again once DM-10153 is finished.
139  old_calibs = []
140  if self.do_photometry:
141  old_calibs = [ref.get('calexp_photoCalib') for ref in data_refs]
142 
143  self.old_dist, self.old_flux, self.old_ref_flux, self.old_source = compute(old_cats, old_calibs)
144 
145  # Update coordinates with the new wcs, and get the new photoCalibs.
146  new_cats = [ref.get('src') for ref in data_refs]
147  new_wcss = []
148  if self.do_astrometry:
149  new_wcss = [ref.get('jointcal_wcs') for ref in data_refs]
150  new_calibs = []
151  if self.do_photometry:
152  new_calibs = [ref.get('jointcal_photoCalib') for ref in data_refs]
153  if self.do_astrometry:
154  for wcs, cat in zip(new_wcss, new_cats):
155  # update in-place the object coordinates based on the new wcs
157 
158  self.new_dist, self.new_flux, self.new_ref_flux, self.new_source = compute(new_cats, new_calibs)
159 
160  if self.verbose:
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))
169 
170  if self.do_photometry:
171  self._photometric_rms()
172  else:
173  self.new_PA1 = None
174 
175  def rms_total(data):
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)
180 
181  if self.do_astrometry:
182  self.old_dist_total = MatchDict(*(tuple(map(rms_total, self.old_dist))*u.radian).to(u.arcsecond))
183  self.new_dist_total = MatchDict(*(tuple(map(rms_total, self.new_dist))*u.radian).to(u.arcsecond))
184  else:
185  self.old_dist_total = MatchDict(None, None)
186  self.new_dist_total = MatchDict(None, None)
187 
188  if self.verbose:
189  if self.do_astrometry:
190  print("relative, absolute astrometry:",
191  self.new_dist_total.relative,
192  self.new_dist_total.absolute)
193  if self.do_photometry:
194  print("Photometric Accuracy (PA1):", self.new_PA1)
195  Rms_result = collections.namedtuple("Rms_result", ["dist_relative", "dist_absolute", "pa1"])
196  return Rms_result(self.new_dist_total.relative, self.new_dist_total.absolute, self.new_PA1)
197 
198  def make_plots(self, data_refs, old_wcs_list,
199  name='', interactive=False, per_ccd_plot=False, outdir='.plots'):
200  """
201  Make plots of various quantites to help with debugging.
202  Requires that `compute_rms()` was run first.
203 
204  Parameters
205  ----------
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.
210  name : str
211  Name to include in plot titles and save files.
212  interactive : bool
213  Turn on matplotlib interactive mode and drop into a debugger when
214  plotting is finished. Otherwise, use a non-interactive backend.
215  per_ccd_plot : bool
216  Plot the WCS per CCD (takes longer and generates many plots for a large camera)
217  outdir : str
218  directory to save plots to.
219  """
220  import matplotlib
221 
222  if not interactive:
223  # Use a non-interactive backend for faster plotting.
224  matplotlib.use('pdf')
225 
226  import matplotlib.pyplot as plt
227  import astropy.visualization
228  # make quantities behave nicely when plotted.
229  astropy.visualization.quantity_support()
230  if interactive:
231  plt.ion()
232 
233  self.log.info("N data_refs: %d", len(data_refs))
234 
235  if self.do_photometry:
236  plot_flux_distributions(plt, self.old_mag, self.new_mag,
238  self.faint, self.bright, self.old_PA1, self.new_PA1,
239  name=name, outdir=outdir)
240  self.log.info("Photometric accuracy (old, new): {:.2e} {:.2e}".format(self.old_PA1,
241  self.new_PA1))
242 
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)
246 
247  if self.do_astrometry:
248  old_dist_rms = MatchDict(*(tuple(map(rms_per_source, self.old_dist))))
249  new_dist_rms = MatchDict(*(tuple(map(rms_per_source, self.new_dist))))
250 
251  self.log.info("relative RMS (old, new): {:.2e} {:.2e}".format(self.old_dist_total.relative,
252  self.new_dist_total.relative))
253  self.log.info("absolute RMS (old, new): {:.2e} {:.2e}".format(self.old_dist_total.absolute,
254  self.new_dist_total.absolute))
255  plot_rms_histogram(plt, old_dist_rms.relative, old_dist_rms.absolute,
256  new_dist_rms.relative, new_dist_rms.absolute,
257  self.old_dist_total.relative, self.old_dist_total.absolute,
258  self.new_dist_total.relative, self.new_dist_total.absolute,
259  name=name, outdir=outdir)
260 
261  plot_all_wcs_deltas(plt, data_refs, self.visits_per_dataRef, old_wcs_list,
262  per_ccd_plot=per_ccd_plot,
263  name=name, outdir=outdir)
264 
265  if interactive:
266  plt.show()
267  import pdb
268  pdb.set_trace()
269 
270  def _photometric_rms(self, sn_cut=300, magnitude_range=3):
271  """
272  Compute the photometric RMS and the photometric repeatablity values (PA1).
273 
274  Parameters
275  ----------
276  sn_cut : float
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.
280  """
281  def rms(flux, ref_flux):
282  return np.sqrt([np.mean((ref_flux[dd] - flux[dd])**2) for dd in flux])
283 
284  self.old_rms = MatchDict(*map(rms, self.old_flux, self.old_ref_flux))
285  self.new_rms = MatchDict(*map(rms, self.new_flux, self.new_ref_flux))
286 
287  # we want to use the absolute fluxes for all of these calculations.
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)
290  self.old_mag = u.Quantity(self.old_ref, u.nJy).to_value(u.ABmag)
291  self.new_mag = u.Quantity(self.new_ref, u.nJy).to_value(u.ABmag)
292 
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])
298  return result
299 
300  old_sn = signal_to_noise(self.old_source.absolute)
301  # Find the faint/bright magnitude limits that are the "flat" part of the rms/magnitude relation.
302  self.faint = self.old_mag[old_sn > sn_cut].max()
303  self.bright = self.faint - magnitude_range
304  if self.verbose:
305  print("PA1 Magnitude range: {:.3f}, {:.3f}".format(self.bright, self.faint))
306  old_good = (self.old_mag < self.faint) & (self.old_mag > self.bright)
307  new_good = (self.new_mag < self.faint) & (self.new_mag > self.bright)
308  self.old_weighted_rms = self.old_rms.absolute/self.old_ref
309  self.new_weighted_rms = self.new_rms.absolute/self.new_ref
310  self.old_PA1 = np.median(self.old_weighted_rms[old_good])
311  self.new_PA1 = np.median(self.new_weighted_rms[new_good])
312 
313  def _make_match_dict(self, reference, visit_catalogs, photoCalibs, refcalib=None):
314  """
315  Return several dicts of sourceID:[values] over the catalogs, to be used in RMS calculations.
316 
317  Parameters
318  ----------
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.
329 
330  Returns
331  -------
332  distances : dict
333  dict of sourceID: array(separation distances for that source)
334  fluxes : dict
335  dict of sourceID: array(fluxes (nJy) for that source)
336  ref_fluxes : dict
337  dict of sourceID: flux (nJy) of the reference object
338  sources : dict
339  dict of sourceID: list(each SourceRecord that was position-matched
340  to this sourceID)
341  """
342  # If we have no photoCalibs, make it the same length as the others for zipping.
343  if photoCalibs == []:
344  photoCalibs = [[]]*len(visit_catalogs)
345 
346  distances = collections.defaultdict(list)
347  fluxes = collections.defaultdict(list)
348  ref_fluxes = {}
349  sources = collections.defaultdict(list)
350  if 'slot_CalibFlux_instFlux' in reference.schema:
351  ref_flux_key = 'slot_CalibFlux'
352  else:
353  ref_flux_key = '{}_flux'
354 
355  def get_fluxes(photoCalib, match):
356  """Return (flux, ref_flux) or None if either is invalid."""
357  # NOTE: Protect against negative fluxes: ignore this match if we find one.
358  flux = match[1]['slot_CalibFlux_instFlux']
359  if flux < 0:
360  return None
361  else:
362  flux = photoCalib.instFluxToNanojansky(match[1], "slot_CalibFlux").value
363 
364  # NOTE: Have to protect against negative "reference" fluxes too.
365  if 'slot' in ref_flux_key:
366  ref_flux = match[0][ref_flux_key+'_instFlux'] # relative reference flux
367  if ref_flux < 0:
368  return None
369  else:
370  ref_flux = refcalib.instFluxToNanojansky(match[0], ref_flux_key).value
371  else:
372  # refcat fluxes are already in nanojansky.
373  ref_flux = match[0][ref_flux_key.format(filt)]
374  if ref_flux < 0:
375  return None
376 
377  Flux = collections.namedtuple('Flux', ('flux', 'ref_flux'))
378  return Flux(flux, ref_flux)
379 
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
382  # things the classifier called sources are not extended.
383  good &= (cat.get('base_ClassificationExtendedness_value') == 0)
384  matches = lsst.afw.table.matchRaDec(reference, cat[good], self.match_radius)
385  for m in matches:
386  if self.do_photometry:
387  flux = get_fluxes(photoCalib, m)
388  if flux is None:
389  continue
390  else:
391  fluxes[m[0].getId()].append(flux.flux)
392  # we can just use assignment here, since the value is always the same.
393  ref_fluxes[m[0].getId()] = flux.ref_flux
394 
395  if self.do_astrometry:
396  # Just use the computed separation distance directly.
397  distances[m[0].getId()].append(m[2])
398 
399  sources[m[0].getId()].append(m[1])
400  # Convert to numpy array for easier math
401  for source in distances:
402  distances[source] = np.array(distances[source])
403  for source in fluxes:
404  fluxes[source] = np.array(fluxes[source])
405 
406  return distances, fluxes, ref_fluxes, sources
407 
408  def _make_visit_catalogs(self, catalogs, visits):
409  """
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.
413 
414  Parameters
415  ----------
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.
420 
421  Returns
422  -------
423  dict
424  dict of visit: catalog of all sources from all CCDs of that visit.
425  """
426  visit_dict = {v: lsst.afw.table.SourceCatalog(catalogs[0].schema) for v in visits}
427  for v, cat in zip(visits, catalogs):
428  visit_dict[v].extend(cat)
429  # We want catalog contiguity to do object selection later.
430  for v in visit_dict:
431  visit_dict[v] = visit_dict[v].copy(deep=True)
432 
433  return visit_dict
434 
435 
436 def plot_flux_distributions(plt, old_mag, new_mag, old_weighted_rms, new_weighted_rms,
437  faint, bright, old_PA1, new_PA1,
438  name='', outdir='.plots'):
439  """Plot various distributions of fluxes and magnitudes.
440 
441  Parameters
442  ----------
443  plt : matplotlib.pyplot instance
444  pyplot instance to plot with
445  old_mag : np.array
446  old magnitudes
447  new_mag : np.array
448  new magnitudes
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))
453  faint : float
454  Faint end of range that PA1 was computed from.
455  bright : float
456  Bright end of range that PA1 was computed from.
457  old_PA1 : float
458  Old value of PA1, to plot as horizontal line.
459  new_PA1 : float
460  New value of PA1, to plot as horizontal line.
461  name : str
462  Name to include in plot titles and save files.
463  outdir : str, optional
464  Directory to write the saved plots to.
465  """
466 
467  import seaborn
468  seaborn.set_style('whitegrid')
469  import scipy.stats
470 
471  old_color = 'blue'
472  new_color = 'red'
473  plt.figure()
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))
486 
487  plt.figure()
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)')
492  plt.ylabel('number')
493  plt.legend(loc='upper right')
494  filename = os.path.join(outdir, '{}-photometry-rms.pdf')
495  plt.savefig(filename.format(name))
496 
497 
498 def plot_all_wcs_deltas(plt, data_refs, visits, old_wcs_list, per_ccd_plot=False,
499  name='', outdir='.plots'):
500  """
501  Various plots of the difference between old and new Wcs.
502 
503  Parameters
504  ----------
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!)
515  name : str
516  Name to include in plot titles and save files.
517  outdir : str, optional
518  Directory to write the saved plots to.
519  """
520 
521  plot_wcs_magnitude(plt, data_refs, visits, old_wcs_list, name, outdir=outdir)
522  plot_all_wcs_quivers(plt, data_refs, visits, old_wcs_list, name, outdir=outdir)
523 
524  if per_ccd_plot:
525  for i, ref in enumerate(data_refs):
526  md = ref.get('calexp_md')
527  dims = lsst.afw.image.bboxFromMetadata(md).getDimensions()
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,
531  outdir=outdir)
532 
533 
534 def make_xy_wcs_grid(x_dim, y_dim, wcs1, wcs2, num=50):
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)
538  x1, y1 = wcs_convert(x, y, wcs1)
539  x2, y2 = wcs_convert(x, y, wcs2)
540  return x1, y1, x2, y2
541 
542 
543 def wcs_convert(xv, yv, wcs):
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()
552  return xout, yout
553 
554 
555 def plot_all_wcs_quivers(plt, data_refs, visits, old_wcs_list, name, outdir='.plots'):
556  """
557  Make quiver plots of the WCS deltas for each CCD in each visit.
558 
559  Parameters
560  ----------
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.
569  name : str
570  Name to include in plot titles and save files.
571  outdir : str, optional
572  Directory to write the saved plots to.
573  """
574 
575  for visit in visits:
576  fig = plt.figure()
577  # fig.set_tight_layout(True)
578  ax = fig.add_subplot(111)
579  for old_wcs, ref in zip(old_wcs_list, data_refs):
580  if ref.dataId['visit'] != visit:
581  continue
582  md = ref.get('calexp_md')
583  dims = lsst.afw.image.bboxFromMetadata(md).getDimensions()
584  Q = plot_wcs_quivers(ax, old_wcs, ref.get('jointcal_wcs'),
585  dims.getX(), dims.getY())
586  # TODO: add CCD bounding boxes to plot once DM-5503 is finished.
587  # TODO: add a circle for the full focal plane.
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')
590  plt.xlabel('RA')
591  plt.ylabel('Dec')
592  plt.title('visit: {}'.format(visit))
593  filename = os.path.join(outdir, '{}-{}-quivers.pdf')
594  plt.savefig(filename.format(name, visit))
595 
596 
597 def plot_wcs_quivers(ax, wcs1, wcs2, x_dim, y_dim):
598  """
599  Plot the delta between wcs1 and wcs2 as vector arrows.
600 
601  Parameters
602  ----------
603  ax : matplotlib.axis
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.
609  x_dim : int
610  Size of array in X-coordinate to make the grid over.
611  y_dim : int
612  Size of array in Y-coordinate to make the grid over.
613  """
614 
615  x1, y1, x2, y2 = make_xy_wcs_grid(x_dim, y_dim, wcs1, wcs2)
616  uu = x2 - x1
617  vv = y2 - y1
618  return ax.quiver(x1, y1, uu, vv, units='x', pivot='tail', scale=1e-3, width=1e-5)
619 
620 
621 def plot_wcs_magnitude(plt, data_refs, visits, old_wcs_list, name, outdir='.plots'):
622  """Plot the magnitude of the WCS change between old and new visits as a heat map.
623 
624  Parameters
625  ----------
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.
634  name : str
635  Name to include in plot titles and save files.
636  outdir : str, optional
637  Directory to write the saved plots to.
638  """
639  for visit in visits:
640  fig = plt.figure()
641  fig.set_tight_layout(True)
642  ax = fig.add_subplot(111)
643  # Start min/max at the "opposite" ends so they always get the first valid value.
644  xmin = np.inf
645  ymin = np.inf
646  xmax = -np.inf
647  ymax = -np.inf
648  for old_wcs, ref in zip(old_wcs_list, data_refs):
649  if ref.dataId['visit'] != visit:
650  continue
651  md = ref.get('calexp_md')
652  dims = lsst.afw.image.bboxFromMetadata(md).getDimensions()
653  x1, y1, x2, y2 = make_xy_wcs_grid(dims.getX(), dims.getY(),
654  old_wcs, ref.get('jointcal_wcs'))
655  uu = x2 - x1
656  vv = y2 - y1
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'))
665  # TODO: add CCD bounding boxes to the plot once DM-5503 is finished.
666  # TODO: add a circle for the full focal plane.
667 
668  # We're reusing only one of the returned images here for colorbar scaling,
669  # but it doesn't matter because we set vmin/vmax so they are all scaled the same.
670  cbar = plt.colorbar(img)
671  cbar.ax.set_ylabel('distortion (arcseconds)')
672  plt.xlim(xmin, xmax)
673  plt.ylim(ymin, ymax)
674  plt.xlabel('RA')
675  plt.ylabel('Dec')
676  plt.title('visit: {}'.format(visit))
677  filename = os.path.join(outdir, '{}-{}-heatmap.pdf')
678  plt.savefig(filename.format(name, visit))
679 
680 
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.
683 
684  Parameters
685  ----------
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.
692  x_dim : int
693  Size of array in X-coordinate to make the grid over.
694  y_dim : int
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.
698  name : str
699  Name to include in plot titles and save files.
700  outdir : str, optional
701  Directory to write the saved plots to.
702  """
703 
704  plt.figure()
705 
706  x1, y1, x2, y2 = make_xy_wcs_grid(x_dim, y_dim, wcs1, wcs2, num=50)
707  plt.plot((x1 - x2) + center[0], (y1 - y2) + center[1], '-')
708  plt.xlabel('delta RA (arcsec)')
709  plt.ylabel('delta Dec (arcsec)')
710  plt.title(name)
711  filename = os.path.join(outdir, '{}-wcs.pdf')
712  plt.savefig(filename.format(name))
713 
714 
715 def plot_rms_histogram(plt, old_rms_relative, old_rms_absolute,
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.
720 
721  Parameters
722  ----------
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
741  name : str
742  Name to include in plot titles and save files.
743  outdir : str, optional
744  Directory to write the saved plots to.
745  """
746  plt.figure()
747 
748  color_rel = 'black'
749  ls_old = 'dotted'
750  color_abs = 'green'
751  ls_new = 'dashed'
752  plotOptions = {'lw': 2, 'range': (0, 0.1)*u.arcsecond, 'normed': True,
753  'bins': 30, 'histtype': 'step'}
754 
755  plt.title('relative vs. absolute: %d vs. %d'%(len(old_rms_relative), len(old_rms_absolute)))
756 
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)
759 
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)
762 
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)
767 
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))
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
lsst::jointcal.utils.plot_all_wcs_quivers
def plot_all_wcs_quivers(plt, data_refs, visits, old_wcs_list, name, outdir='.plots')
Definition: utils.py:555
lsst::log.log.logContinued.info
def info(fmt, *args)
Definition: logContinued.py:198
lsst::jointcal.utils.JointcalStatistics.old_dist_total
old_dist_total
Definition: utils.py:182
lsst::jointcal.utils.JointcalStatistics.compute_rms
def compute_rms(self, data_refs, reference)
Definition: utils.py:85
lsst::jointcal.utils.MatchDict
MatchDict
Definition: utils.py:46
lsst::jointcal.utils.JointcalStatistics.new_source
new_source
Definition: utils.py:158
lsst::jointcal.utils.JointcalStatistics.do_astrometry
do_astrometry
Definition: utils.py:79
lsst::jointcal.utils.JointcalStatistics.log
log
Definition: utils.py:81
pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
ast::append
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
lsst::jointcal.utils.plot_wcs
def plot_wcs(plt, wcs1, wcs2, x_dim, y_dim, center=(0, 0), name="", outdir='.plots')
Definition: utils.py:681
lsst::jointcal.utils.JointcalStatistics.do_photometry
do_photometry
Definition: utils.py:78
lsst::log::Log::getLogger
static Log getLogger(Log const &logger)
Definition: Log.h:760
std::tuple
lsst::jointcal.utils.plot_all_wcs_deltas
def plot_all_wcs_deltas(plt, data_refs, visits, old_wcs_list, per_ccd_plot=False, name='', outdir='.plots')
Definition: utils.py:498
lsst::jointcal.utils.JointcalStatistics.old_mag
old_mag
Definition: utils.py:290
lsst::jointcal.utils.JointcalStatistics.new_weighted_rms
new_weighted_rms
Definition: utils.py:309
lsst::jointcal.utils.JointcalStatistics.match_radius
match_radius
Definition: utils.py:76
lsst::jointcal.utils.JointcalStatistics.verbose
verbose
Definition: utils.py:80
lsst::jointcal.utils.JointcalStatistics.bright
bright
Definition: utils.py:303
lsst::jointcal.utils.JointcalStatistics._photometric_rms
def _photometric_rms(self, sn_cut=300, magnitude_range=3)
Definition: utils.py:270
lsst::jointcal.utils.JointcalStatistics.make_plots
def make_plots(self, data_refs, old_wcs_list, name='', interactive=False, per_ccd_plot=False, outdir='.plots')
Definition: utils.py:198
lsst::jointcal.utils.JointcalStatistics.filters
filters
Definition: utils.py:108
lsst::afw::table::updateSourceCoords
void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList)
Update sky coordinates in a collection of source objects.
Definition: wcsUtils.cc:96
lsst::jointcal.utils.JointcalStatistics.old_source
old_source
Definition: utils.py:143
lsst::afw::table._source.SourceCatalog
Definition: _source.py:33
lsst::log
Definition: Log.h:706
lsst::jointcal.utils.JointcalStatistics.new_dist_total
new_dist_total
Definition: utils.py:183
lsst::afw::table::matchRaDec
template SourceMatchVector matchRaDec(SourceCatalog const &, lsst::geom::Angle, MatchControl const &)
max
int max
Definition: BoundedField.cc:104
lsst::afw::table
Definition: table.dox:3
std::map
STL class.
lsst::jointcal.utils.plot_wcs_magnitude
def plot_wcs_magnitude(plt, data_refs, visits, old_wcs_list, name, outdir='.plots')
Definition: utils.py:621
lsst::jointcal.utils.JointcalStatistics.old_weighted_rms
old_weighted_rms
Definition: utils.py:308
lsst::jointcal.utils.plot_wcs_quivers
def plot_wcs_quivers(ax, wcs1, wcs2, x_dim, y_dim)
Definition: utils.py:597
lsst::jointcal.utils.JointcalStatistics.new_ref
new_ref
Definition: utils.py:289
to
table::Key< int > to
Definition: TransformMap.cc:349
lsst::meas::base::Flux
double Flux
Definition: constants.h:51
lsst::geom
Definition: geomOperators.dox:4
lsst::afw::image::bboxFromMetadata
lsst::geom::Box2I bboxFromMetadata(daf::base::PropertySet &metadata)
Determine the image bounding box from its metadata (FITS header)
Definition: Image.cc:694
lsst::jointcal.utils.JointcalStatistics._make_match_dict
def _make_match_dict(self, reference, visit_catalogs, photoCalibs, refcalib=None)
Definition: utils.py:313
lsst::jointcal.utils.JointcalStatistics.new_rms
new_rms
Definition: utils.py:285
lsst::jointcal.utils.JointcalStatistics.old_PA1
old_PA1
Definition: utils.py:310
min
int min
Definition: BoundedField.cc:103
lsst::jointcal.utils.JointcalStatistics.__init__
def __init__(self, match_radius=0.1 *arcseconds, flux_limit=100.0, do_photometry=True, do_astrometry=True, verbose=False)
Definition: utils.py:60
lsst::jointcal.utils.JointcalStatistics.old_rms
old_rms
Definition: utils.py:284
lsst::jointcal.utils.plot_flux_distributions
def plot_flux_distributions(plt, old_mag, new_mag, old_weighted_rms, new_weighted_rms, faint, bright, old_PA1, new_PA1, name='', outdir='.plots')
Definition: utils.py:436
lsst::jointcal.utils.make_xy_wcs_grid
def make_xy_wcs_grid(x_dim, y_dim, wcs1, wcs2, num=50)
Definition: utils.py:534
lsst::jointcal.utils.JointcalStatistics.old_ref
old_ref
Definition: utils.py:288
lsst::jointcal.utils.JointcalStatistics.visits_per_dataRef
visits_per_dataRef
Definition: utils.py:109
lsst::jointcal.utils.JointcalStatistics.new_mag
new_mag
Definition: utils.py:291
lsst::jointcal.utils.JointcalStatistics.new_PA1
new_PA1
Definition: utils.py:173
lsst::jointcal.utils.JointcalStatistics._make_visit_catalogs
def _make_visit_catalogs(self, catalogs, visits)
Definition: utils.py:408
lsst::jointcal.utils.JointcalStatistics.faint
faint
Definition: utils.py:302
lsst::jointcal.utils.JointcalStatistics
Definition: utils.py:49
lsst::jointcal.utils.JointcalStatistics.flux_limit
flux_limit
Definition: utils.py:77
lsst::jointcal.utils.plot_rms_histogram
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')
Definition: utils.py:715
lsst::jointcal.utils.wcs_convert
def wcs_convert(xv, yv, wcs)
Definition: utils.py:543