LSSTApplications  20.0.0
LSSTDataManagementBasePackage
matchFakes.py
Go to the documentation of this file.
1 """
2 matchFakes.py.
3 
4 Matches fakes based on position stored in the calibrated exposure image header
5 """
6 
7 import re
8 import argparse
9 import collections
10 
11 import numpy as np
12 import astropy.table
13 
14 import lsst.daf.persistence as dafPersist
15 import lsst.afw.geom.ellipses
16 import lsst.afw.geom
17 import lsst.pex.exceptions
18 from lsst.afw.table import SourceCatalog, SchemaMapper
19 
20 NO_FOOTPRINT = lsst.afw.table.SOURCE_IO_NO_FOOTPRINTS
21 
22 
23 def combineWithForce(meas, force):
24  """Combine the meas and forced_src catalogs."""
25  if len(meas) != len(force):
26  raise Exception("# Meas and Forced_src catalogs should have "
27  "the same size!")
28  mapper = SchemaMapper(meas.schema)
29  mapper.addMinimalSchema(meas.schema)
30  newSchema = mapper.getOutputSchema()
31  # Add new fields
32  newSchema.addField('force_deblend_nChild', type=np.int32)
33  newSchema.addField('force_base_ClassificationExtendedness_value', type=float)
34  newSchema.addField('force_ext_photometryKron_KronFlux_instFlux', type=float)
35  newSchema.addField('force_ext_photometryKron_KronFlux_instFluxErr', type=float)
36  newSchema.addField('force_base_PsfFlux_instFlux', type=float)
37  newSchema.addField('force_base_PsfFlux_instFluxErr', type=float)
38  newSchema.addField('force_ext_photometryKron_KronFlux_apCorr', type=float)
39  newSchema.addField('force_ext_photometryKron_KronFlux_apCorrErr', type=float)
40  newSchema.addField('force_base_PsfFlux_apCorr', type=float)
41  newSchema.addField('force_base_PsfFlux_apCorrErr', type=float)
42  newSchema.addField('force_modelfit_CModel_instFlux', type=float)
43  newSchema.addField('force_modelfit_CModel_instFluxErr', type=float)
44  newSchema.addField('force_modelfit_CModel_fracDev', type=float)
45  newSchema.addField('force_modelfit_CModel_exp_instFlux', type=float)
46  newSchema.addField('force_modelfit_CModel_exp_instFluxErr', type=float)
47  newSchema.addField('force_modelfit_CModel_dev_instFlux', type=float)
48  newSchema.addField('force_modelfit_CModel_dev_instFluxErr', type=float)
49  newSchema.addField('force_modelfit_CModel_apCorr', type=float)
50  newSchema.addField('force_modelfit_CModel_apCorrErr', type=float)
51  newSchema.addField('force_modelfit_CModel_exp_apCorr', type=float)
52  newSchema.addField('force_modelfit_CModel_exp_apCorrErr', type=float)
53  newSchema.addField('force_modelfit_CModel_dev_apCorr', type=float)
54  newSchema.addField('force_modelfit_CModel_dev_apCorrErr', type=float)
55 
56  newCols = ['deblend_nChild', 'base_ClassificationExtendedness_value',
57  'ext_photometryKron_KronFlux_instFlux', 'ext_photometryKron_KronFlux_instFluxErr',
58  'base_PsfFlux_instFlux', 'base_PsfFlux_instFluxErr',
59  'ext_photometryKron_KronFlux_apCorr', 'ext_photometryKron_KronFlux_apCorrErr',
60  'base_PsfFlux_apCorr', 'base_PsfFlux_apCorrErr',
61  'modelfit_CModel_instFlux', 'modelfit_CModel_instFluxErr',
62  'modelfit_CModel_exp_apCorr', 'modelfit_CModel_exp_apCorrErr',
63  'modelfit_CModel_exp_instFlux', 'modelfit_CModel_exp_instFlux',
64  'modelfit_CModel_exp_apCorr', 'modelfit_CModel_exp_apCorrErr',
65  'modelfit_CModel_dev_instFlux', 'modelfit_CModel_dev_instFluxErr',
66  'modelfit_CModel_dev_apCorr', 'modelfit_CModel_dev_apCorrErr',
67  'modelfit_CModel_fracDev']
68  measAlias = meas.schema.getAliasMap()
69  newAlias = newSchema.getAliasMap()
70  for aliasKey in measAlias.keys():
71  newAlias.set(aliasKey, measAlias[aliasKey])
72  combSrc = SourceCatalog(newSchema)
73  combSrc.extend(meas, mapper=mapper)
74 
75  for key in newCols:
76  combSrc['force_' + key][:] = force[key][:]
77 
78  return combSrc
79 
80 
81 def getMag(flux, fluxerr, zeropoint):
82  """Return the magnitude and error."""
83  mag = -2.5 * np.log10(flux)
84  magerr = (2.5 / np.log(10.0) * fluxerr / flux)
85  return (mag.T + zeropoint).T, magerr
86 
87 
88 def getEllipse(quad):
89  """Return the Re, b/a and PA for a given quadrupole moment."""
90  e = lsst.afw.geom.ellipses.Axes(quad)
91  rad = e.getA()
92  try:
93  b_a = (e.getB() / e.getA())
94  except ZeroDivisionError:
95  b_a = np.nan
96  pa = (e.getTheta() * 180.0 / np.pi)
97  return rad, b_a, pa
98 
99 
100 def matchToFakeCatalog(sources, fakeCatalog):
101  """
102  Match to the fake catalog and append those columns to the source table.
103 
104  this assumes the sources are an astropy table
105  or it will throw a TypeError
106  """
107  if not isinstance(sources, astropy.table.Table):
108  raise TypeError("Expect an astropy table for sources"
109  " use getAstroTable to convert")
110 
111  fakes = astropy.table.Table().read(fakeCatalog)
112  fakes.rename_column('ID', 'fakeId')
113  return astropy.table.join(sources, fakes, keys='fakeId', join_type='left')
114 
115 
116 def getFakeMatchesHeader(cal_md, sources, tol=1.0):
117  """
118  Return the fake matches based on the information in the header.
119 
120  returns a tuple with:
121  the positions in pixels of the fake sources added to the chip
122  the match is in a dictionary of the form:
123  {fakeid:[ind_of_match_in_sources,...],...}
124 
125  look within a tolerance of 1 pixel in each direction
126  """
127  fakeXY = collections.defaultdict(tuple)
128  fakename = re.compile('FAKE([0-9]+)')
129  for card in cal_md.names():
130  m = fakename.match(card)
131  if m is not None:
132  x, y = list(map(float, (cal_md.getScalar(card)).split(',')))
133  fakeXY[int(m.group(1))] = (x, y)
134 
135  srcX, srcY = sources.getX(), sources.getY()
136  srcIndex = collections.defaultdict(list)
137  for fid, fcoord in fakeXY.items():
138  distX = srcX - fcoord[0]
139  distY = srcY - fcoord[1]
140  matched = (np.abs(distX) < tol) & (np.abs(distY) < tol)
141  srcIndex[fid] = np.where(matched)[0]
142 
143  return fakeXY, srcIndex
144 
145 
146 def getFakeMatchesRaDec(sources, radecCatFile, bbox, wcs, tol=1.0,
147  reffMatch=False, pix=0.168, minRad=None,
148  raCol='RA', decCol='Dec'):
149  """
150  Return the fake matches based on an radec match.
151 
152  Args:
153  sources: source table object
154  radecCatFile: filename for fits file of fake object table,
155  including ra/dec
156  bbox: Bounding Box of exposure (ccd or patch) within which
157  to match fakes
158  wcs: Wcs of source image
159 
160  KeywordArgs:
161  tol: tolerance within which to match sources, given in PIXELS
162 
163  Returns:
164  fakeXY: set of pixel positions of fake sources
165  srcIdx: matches to fake sources in a dictionary of the form
166  {fakeid:[ind_of_match_in_sources,...],...}
167 
168  Raise:
169  IOError: couldn't open radecCatFile
170  """
171  fakeXY = collections.defaultdict(tuple)
172  try:
173  fakeCat = astropy.table.Table().read(radecCatFile, format='fits')
174  except IOError:
175  raise
176 
177  for fakeSrc in fakeCat:
178 
179  fakePix = lsst.afw.geom.SpherePoint(fakeSrc[raCol], fakeSrc[decCol], lsst.afw.geom.degrees)
180  fakeCoord = wcs.skyToPixel(fakePix)
181 
182  if bbox.contains(fakeCoord):
183  if reffMatch:
184  fakeXY[int(fakeSrc['ID'])] = (fakeCoord.getX(),
185  fakeCoord.getY(),
186  (fakeSrc['reff'] / pix) * tol)
187  else:
188  fakeXY[int(fakeSrc['ID'])] = (fakeCoord.getX(),
189  fakeCoord.getY(),
190  tol)
191 
192  srcX, srcY = sources.getX(), sources.getY()
193  srcIndex = collections.defaultdict(list)
194  srcClose = collections.defaultdict(list)
195  for fid, fcoord in fakeXY.items():
196  distX = (srcX - fcoord[0])
197  distY = (srcY - fcoord[1])
198  distR = np.sqrt((np.abs(distX) ** 2.0) + (np.abs(distY) ** 2.0))
199  closest = np.nanargmin(distR)
200  radMatch = fcoord[2]
201  if minRad is not None:
202  if radMatch < minRad:
203  radMatch = minRad
204  if reffMatch:
205  matched = (distR <= radMatch)
206  else:
207  matched = (np.abs(distX) <= radMatch) & (np.abs(distY) <= radMatch)
208 
209  srcIndex[fid] = np.where(matched)[0]
210  srcClose[fid] = closest
211 
212  return fakeXY, srcIndex, srcClose
213 
214 
215 def getFakeSources(butler, dataId, tol=1.0,
216  extraCols=('zeropoint', 'visit', 'ccd'),
217  includeMissing=False, footprints=False, radecMatch=None,
218  multiband=False, reffMatch=False, pix=0.168, minRad=None,
219  raCol='RA', decCol='Dec'):
220  """
221  Get list of sources which agree in pixel position with fake ones with tol.
222 
223  This returns a sourceCatalog of all the matched fake objects,
224  note, there will be duplicates in this list, since I haven't
225  checked deblend.nchild, and I'm only doing a tolerance match,
226  which could include extra sources
227 
228  The outputs can include extraCols as long as they are one of:
229  zeropoint, visit, ccd, thetaNorth, pixelScale
230 
231  If includeMissing is true, then the pipeline looks at the fake sources
232  added in the header and includes an entry in the table for sources without
233  any measurements, specifically the 'id' column will be 0
234 
235  radecMatch is the fakes table. if it's not None(default), then do an ra/dec
236  match with the input catalog instead of looking in the header for where the
237  sources where added
238  """
239  coaddData = "deepCoadd_calexp"
240  coaddMeta = "deepCoadd_calexp_md"
241 
242  availExtras = {'zeropoint': {'type': float, 'doc': 'zeropoint'},
243  'visit': {'type': int, 'doc': 'visit id'},
244  'ccd': {'type': int, 'doc': 'ccd id'},
245  'thetaNorth': {'type': lsst.afw.geom.Angle,
246  'doc': 'angle to north'},
247  'pixelScale': {'type': float,
248  'doc': 'pixelscale in arcsec/pixel'}}
249 
250  if not np.in1d(extraCols, list(availExtras.keys())).all():
251  print("extraCols must be in ", availExtras)
252 
253  try:
254  if 'filter' not in dataId:
255  sources = butler.get('src', dataId,
256  flags=lsst.afw.table.SOURCE_IO_NO_FOOTPRINTS,
257  immediate=True)
258  cal = butler.get('calexp', dataId, immediate=True)
259  cal_md = butler.get('calexp_md', dataId, immediate=True)
260  else:
261  meas = butler.get('deepCoadd_meas', dataId,
262  flags=NO_FOOTPRINT,
263  immediate=True)
264  force = butler.get('deepCoadd_forced_src', dataId,
265  flags=NO_FOOTPRINT,
266  immediate=True)
267  sources = combineWithForce(meas, force)
268  cal = butler.get(coaddData, dataId, immediate=True)
269  cal_md = butler.get(coaddMeta, dataId, immediate=True)
270  except RuntimeError:
271  print("skipping", dataId)
272  return None
273 
274  if ('pixelScale' in extraCols) or ('thetaNorth' in extraCols):
275  wcs = cal.getWcs()
276  availExtras['pixelScale']['value'] = wcs.getPixelScale().asArcseconds()
277  # The 8 lines of code below find the angle to north, first the mid pixel of the calexp is found,
278  # then the pixel to sky matrix at this point, the coordinate this gives can then be used to find the
279  # linearized sky to pixel matrix which can then be used to find the angle.
280  xMid = cal.getWidth() // 2
281  yMid = cal.getHeight() // 2
282  midPoint = lsst.afw.geom.Point2D(xMid, yMid)
283  midCoord = wcs.pixelToSky(midPoint)
284  northSkyToPixelMatrix = wcs.linearizeSkyToPixel(midCoord, lsst.afw.geom.degrees)
285  northSkyToPixelMatrix = northSkyToPixelMatrix.getLinear()
286  availExtras['thetaNorth']['value'] = (np.arctan2(*tuple(northSkyToPixelMatrix
287  (lsst.afw.geom.Point2D(1.0, 0.0)))))*lsst.afw.geom.radians
288 
289  if 'visit' in extraCols:
290  availExtras['visit']['value'] = dataId['visit']
291  if 'ccd' in extraCols:
292  availExtras['ccd']['value'] = dataId['ccd']
293  if 'zeropoint' in extraCols:
294  zeropoint = 2.5 * np.log10(cal_md.getScalar('FLUXMAG0'))
295  availExtras['zeropoint']['value'] = zeropoint
296 
297  if radecMatch is None:
298  fakeXY, srcIndex = getFakeMatchesHeader(cal_md, sources, tol=tol)
299  else:
300  if minRad is not None:
301  print("# The min matching radius is %4.1f pixel" % minRad)
302  bbox = lsst.afw.geom.Box2D(cal.getBBox(lsst.afw.image.PARENT))
303  fakeXY, srcIndex, srcClose = getFakeMatchesRaDec(sources,
304  radecMatch,
305  bbox,
306  cal.getWcs(),
307  tol=tol,
308  reffMatch=reffMatch,
309  pix=pix,
310  minRad=minRad,
311  raCol=raCol,
312  decCol=decCol)
313 
314  mapper = SchemaMapper(sources.schema)
315  mapper.addMinimalSchema(sources.schema)
316  newSchema = mapper.getOutputSchema()
317  newSchema.addField('fakeId', type=np.int32,
318  doc='id of fake source matched to position')
319  newSchema.addField('nMatched', type=np.int32,
320  doc='Number of matched objects')
321  newSchema.addField('nPrimary', type=np.int32,
322  doc='Number of unique matched objects')
323  newSchema.addField('nNoChild', type=np.int32,
324  doc='Number of matched objects with nchild==0')
325  newSchema.addField('rMatched', type=float,
326  doc='Radius used form atching obects, in pixel')
327  newSchema.addField('fakeOffX', type=float,
328  doc='offset from input fake position in X (pixels)')
329  newSchema.addField('fakeOffY', type=float,
330  doc='offset from input fake position in Y (pixels)')
331  newSchema.addField('fakeOffR', type=float,
332  doc='offset from input fake position in radius')
333  newSchema.addField('fakeClosest', type="Flag",
334  doc='Is this match the closest one?')
335 
336  for extraName in set(extraCols).intersection(availExtras):
337  newSchema.addField(extraName, type=availExtras[extraName]['type'],
338  doc=availExtras[extraName]['doc'])
339 
340  srcList = SourceCatalog(newSchema)
341  srcList.reserve(sum([len(s) for s in srcIndex.values()]) +
342  (0 if not includeMissing else list(srcIndex.values()).count([])))
343 
344  centroidKey = sources.getCentroidKey()
345  isPrimary = sources.schema.find('detect_isPrimary').getKey()
346  nChild = sources.schema.find('force_deblend_nChild').getKey()
347  for ident, sindlist in srcIndex.items():
348  rMatched = fakeXY[ident][2]
349  if minRad is not None:
350  if rMatched < minRad:
351  rMatched = minRad
352  nMatched = len(sindlist)
353  nPrimary = np.sum([sources[int(obj)].get(isPrimary) for obj in sindlist])
354  nNoChild = np.sum([(sources[int(obj)].get(nChild) == 0) for
355  obj in sindlist])
356  if includeMissing and (nMatched == 0):
357  newRec = srcList.addNew()
358  newRec.set('fakeId', ident)
359  newRec.set('id', 0)
360  newRec.set('nMatched', 0)
361  newRec.set('rMatched', rMatched)
362  for ss in sindlist:
363  newRec = srcList.addNew()
364  newRec.assign(sources[int(ss)], mapper)
365  newRec.set('fakeId', ident)
366  newRec.set('nMatched', nMatched)
367  newRec.set('nPrimary', nPrimary)
368  newRec.set('nNoChild', nNoChild)
369  newRec.set('rMatched', rMatched)
370  offsetX = (sources[int(ss)].get(centroidKey).getX() -
371  fakeXY[ident][0])
372  newRec.set('fakeOffX', offsetX)
373  offsetY = (sources[int(ss)].get(centroidKey).getY() -
374  fakeXY[ident][1])
375  newRec.set('fakeOffY', offsetY)
376  newRec.set('fakeOffR', np.sqrt(offsetX ** 2.0 + offsetY ** 2.0))
377  if radecMatch:
378  if int(ss) == int(srcClose[ident]):
379  newRec.set('fakeClosest', True)
380  else:
381  newRec.set('fakeClosest', False)
382 
383  if includeMissing:
384  srcList = srcList.copy(deep=True)
385 
386  for extraName in set(extraCols).intersection(availExtras):
387  tempCol = srcList.get(extraName)
388  tempCol.fill(availExtras[extraName]['value'])
389 
390  return srcList
391 
392 
393 def getAstroTable(src, mags=True):
394  """
395  Return an astropy table with all the src entries.
396 
397  if the entries are complex objects, it breaks them down:
398  ellipse entries are broken into
399  ellipse_a = semi-major axis
400  ellipse_q = axis ratio (always < 1)
401  ellipse_theta = rotation of semi-major axis
402  from chip x-axis in degrees
403  if mags is True, returns the magnitudes for all the flux columns
404  """
405  tab = astropy.table.Table()
406  for name in src.schema.getNames():
407  # For reasons I don't understand a lookup by name is much
408  # slower than a lookup by key
409  nameKey = src.schema.find(name).getKey()
410  try:
411  tab.add_column(astropy.table.Column(name=name,
412  data=src.get(nameKey)))
413  except lsst.pex.exceptions.LsstException:
414  quadrupole = lsst.afw.geom.ellipses.ellipsesLib.Quadrupole
415  if type(src[0].get(nameKey)) is quadrupole:
416  """Check for shape measurements"""
417  reff, q, theta = list(zip(*[getEllipse(s.get(nameKey))
418  for s in src]))
419  tab.add_column(astropy.table.Column(name=name+'_a',
420  data=reff))
421  tab.add_column(astropy.table.Column(name=name+'_q',
422  data=q))
423  tab.add_column(astropy.table.Column(name=name+'_theta',
424  data=theta))
425  elif type(src[0].get(nameKey)) is lsst.afw.geom.SpherePoint:
426  """Check for coordinate measurements"""
427  x, y = list(zip(*[(s.get(nameKey).getRa().asDegrees(),
428  s.get(nameKey).getDec().asDegrees()) for s in src]))
429  tab.add_column(astropy.table.Column(name=name+'_ra', data=x))
430  tab.add_column(astropy.table.Column(name=name+'_dec', data=y))
431  else:
432  keyData = np.array([s.get(nameKey) for s in src])
433  tab.add_column(astropy.table.Column(name=name,
434  data=keyData))
435 
436  if isinstance(src[0].get(nameKey), lsst.afw.geom.Angle):
437  # Report angles in degrees
438  tab.remove_column(name)
439  newCol = astropy.table.Column(data=[s.get(nameKey).asDegrees()
440  for s in src],
441  dtype=float, name=name)
442  tab.add_column(newCol)
443 
444  if mags:
445  # This is a horrible hack, but I don't think we can use the slots,
446  # since not all the fluxes end up in the slots
447  for col in tab.colnames:
448  colMatch = (re.match(r'^flux\.[a-z]+$', col) or
449  re.match(r'^flux\.[a-z]+.apcorr$', col) or
450  re.match(r'^force.flux\.[a-z]+$', col) or
451  re.match(r'^force.flux\.[a-z]+.apcorr$', col) or
452  re.match('^force.cmodel.+flux$', col) or
453  re.match('^force.cmodel.+flux.apcorr$', col) or
454  re.match('^cmodel.+flux$', col) or
455  re.match('^cmodel.+flux.apcorr$', col))
456  if colMatch:
457  zp = tab['zeropoint'] if not re.search('apcorr', col) else 0.0
458  mag, magerr = getMag(tab[col], tab[col+'.err'], zp)
459  tab.add_column(astropy.table.Column(name=re.sub('instFlux', 'mag',
460  col),
461  data=mag))
462  tab.add_column(astropy.table.Column(name=re.sub('instFlux', 'mag',
463  col+'.err'),
464  data=magerr))
465 
466  return tab
467 
468 
469 def returnMatchSingle(butler, slist, visit, ccd,
470  filt=None, tol=1.0, pix=0.168,
471  fakeCat=None, pixMatch=False, multiband=False,
472  reffMatch=False, includeMissing=True, minRad=None,
473  raCol='RA', decCol='Dec'):
474  """Return matched catalog for each CCD or Patch."""
475  if filt is None:
476  print('Doing ccd %d' % int(ccd))
477  mlis = getFakeSources(butler,
478  {'visit': visit, 'ccd': int(ccd)},
479  includeMissing=includeMissing,
480  extraCols=('visit', 'ccd',
481  'zeropoint', 'pixelScale',
482  'thetaNorth'),
483  radecMatch=fakeCat if not pixMatch else None,
484  tol=tol, reffMatch=reffMatch, pix=pix,
485  minRad=minRad, raCol=raCol, decCol=decCol)
486  else:
487  print('Doing patch %s' % ccd)
488  mlis = getFakeSources(butler,
489  {'tract': visit, 'patch': ccd,
490  'filter': filt},
491  includeMissing=includeMissing,
492  extraCols=('thetaNorth', 'pixelScale',
493  'zeropoint'),
494  radecMatch=fakeCat if not pixMatch else None,
495  tol=tol, multiband=multiband,
496  reffMatch=reffMatch, pix=pix,
497  minRad=minRad, raCol=raCol, decCol=decCol)
498 
499  if mlis is None:
500  print(' No match returns!')
501  else:
502  if slist is None:
503  slist = mlis.copy(True)
504  else:
505  slist.extend(mlis, True)
506  del mlis
507 
508  return slist
509 
510 
511 def returnMatchTable(rootDir, visit, ccdList, outfile=None, fakeCat=None,
512  overwrite=False, filt=None, tol=1.0, pixMatch=False,
513  multiband=False, reffMatch=False, pix=0.168,
514  multijobs=1, includeMissing=True, minRad=None,
515  raCol='RA', decCol='Dec'):
516  """
517  Driver (main function) for return match to fakes.
518 
519  INPUT: rootDir = rerun directory
520  visit = visit id (int) (or tracts)
521  ccdList = list of ccds to look at (or patches)
522  outdir = output directory for matched file,
523  None means no output written
524  fakeCat = fake catalog to match to,
525  None means the fake sources are just
526  extracted from the header of the CCDs based on
527  position but no matching is done
528  overwrite = whether to overwrite the existing output file,
529  default is False
530  pixMatch = do pixel matching instead of ra/dec matching
531  even if there is a catalog supplied
532  multiband = whether match to forced photometry catalogs
533  from multiband process
534  reffMatch = whether match fake sources in pixel radius
535  or using tol x Reff (Only for Ra, Dec match)
536  OUTPUT: returns an astropy.table.Table with all the entries
537  from the source catalog for objects which match in pixel
538  position to the fake sources
539  """
540  butler = dafPersist.Butler(rootDir)
541  slist = None
542 
543  if multijobs > 1:
544  try:
545  from joblib import Parallel, delayed
546  mlist = Parallel(n_jobs=multijobs)(delayed(returnMatchSingle)(
547  butler, None, visit, ccd,
548  filt=filt,
549  fakeCat=fakeCat,
550  includeMissing=includeMissing,
551  pixMatch=pixMatch,
552  reffMatch=reffMatch, tol=tol,
553  multiband=multiband,
554  minRad=minRad,
555  pix=pix,
556  decCol=decCol,
557  raCol=raCol) for ccd in ccdList)
558  for m in mlist:
559  if m is not None:
560  if slist is None:
561  slist = m.copy(True)
562  else:
563  slist.extend(m, True)
564  del m
565  except ImportError:
566  print("# Can not import joblib, stop multiprocessing!")
567  for ccd in ccdList:
568  slist = returnMatchSingle(butler, slist, visit, ccd,
569  filt=filt, fakeCat=fakeCat,
570  includeMissing=includeMissing,
571  pixMatch=pixMatch,
572  reffMatch=reffMatch,
573  tol=tol, pix=pix,
574  multiband=multiband,
575  minRad=minRad,
576  raCol=raCol, decCol=decCol)
577  else:
578  for ccd in ccdList:
579  slist = returnMatchSingle(butler, slist, visit, ccd,
580  filt=filt, fakeCat=fakeCat,
581  includeMissing=includeMissing,
582  pixMatch=pixMatch,
583  reffMatch=reffMatch,
584  tol=tol, pix=pix,
585  multiband=multiband,
586  minRad=minRad,
587  raCol=raCol, decCol=decCol)
588 
589  if slist is None:
590  print("Returns no match....!")
591 
592  return None
593  else:
594  astroTable = getAstroTable(slist, mags=True)
595 
596  if fakeCat is not None:
597  astroTable = matchToFakeCatalog(astroTable, fakeCat)
598 
599  if outfile is not None:
600  try:
601  astroTable.write(outfile+'.fits', format='fits',
602  overwrite=overwrite)
603  except IOError:
604  print("Try setting the option -w to overwrite the file.")
605  raise
606 
607  return astroTable
608 
609 
610 if __name__ == '__main__':
611  # TODO: this should use the LSST/HSC conventions
612  parser = argparse.ArgumentParser()
613  parser.add_argument('rootDir', help='root dir of data repo')
614  parser.add_argument('visit',
615  help='id of visit (or tract, if filter is specified)',
616  type=int)
617  parser.add_argument('-f', '--filter', dest='filt',
618  help='name of filter, if none assume single visit',
619  default=None)
620  parser.add_argument('--ccd', nargs='+', help='id of ccd(s) or patches')
621  parser.add_argument('-o', help='outputfilename', default=None,
622  dest='outfile')
623  parser.add_argument('-c', help='fake catalog', default=None,
624  dest='fakeCat')
625  parser.add_argument('-w', '--overwrite', help='overwrite output file',
626  dest='ow', default=False, action='store_true')
627  parser.add_argument('-m', '--multiband',
628  help='Match multiband measurements',
629  dest='multiband', default=False, action='store_true')
630  parser.add_argument('-r', '--reffMatch',
631  help='Match the fake sources using tol x Reff',
632  dest='reffMatch', default=False, action='store_true')
633  parser.add_argument('--min', '--minRad',
634  help='Minimum matching radius in unit of pixel',
635  dest='minRad', type=float, default=None)
636  parser.add_argument('-t', '--tolerance', type=float, dest='tol',
637  help='matching radius in PIXELS (default=1.0)')
638  parser.add_argument('-j', '--multijobs', type=int,
639  help='Number of jobs run at the same time',
640  dest='multijobs', default=1)
641  parser.add_argument('--ra', '--raCol', dest='raCol',
642  help='Name of the column for RA',
643  default='RA')
644  parser.add_argument('--dec', '--decCol', dest='decCol',
645  help='Name of the column for Dec',
646  default='Dec')
647  args = parser.parse_args()
648 
649  returnMatchTable(args.rootDir, args.visit, args.ccd, args.outfile,
650  args.fakeCat, overwrite=args.ow, filt=args.filt,
651  tol=args.tol, multiband=args.multiband,
652  reffMatch=args.reffMatch,
653  multijobs=args.multijobs,
654  minRad=args.minRad,
655  raCol=args.raCol, decCol=args.decCol)
lsst.synpipe.matchFakes.returnMatchSingle
def returnMatchSingle(butler, slist, visit, ccd, filt=None, tol=1.0, pix=0.168, fakeCat=None, pixMatch=False, multiband=False, reffMatch=False, includeMissing=True, minRad=None, raCol='RA', decCol='Dec')
Definition: matchFakes.py:469
lsst.synpipe.matchFakes.combineWithForce
def combineWithForce(meas, force)
Definition: matchFakes.py:23
lsst::afw::geom::ellipses.quadrupole.quadrupoleContinued.Quadrupole
Definition: quadrupoleContinued.py:8
lsst.synpipe.matchFakes.getFakeSources
def getFakeSources(butler, dataId, tol=1.0, extraCols=('zeropoint', 'visit', 'ccd'), includeMissing=False, footprints=False, radecMatch=None, multiband=False, reffMatch=False, pix=0.168, minRad=None, raCol='RA', decCol='Dec')
Definition: matchFakes.py:215
lsst.synpipe.matchFakes.matchToFakeCatalog
def matchToFakeCatalog(sources, fakeCatalog)
Definition: matchFakes.py:100
lsst.synpipe.matchFakes.getFakeMatchesHeader
def getFakeMatchesHeader(cal_md, sources, tol=1.0)
Definition: matchFakes.py:116
lsst::afw::geom::ellipses::Axes
An ellipse core for the semimajor/semiminor axis and position angle parametrization (a,...
Definition: Axes.h:47
lsst.synpipe.matchFakes.getAstroTable
def getAstroTable(src, mags=True)
Definition: matchFakes.py:393
lsst::geom::all
bool all(CoordinateExpr< N > const &expr) noexcept
Return true if all elements are true.
Definition: CoordinateExpr.h:81
lsst.synpipe.matchFakes.getFakeMatchesRaDec
def getFakeMatchesRaDec(sources, radecCatFile, bbox, wcs, tol=1.0, reffMatch=False, pix=0.168, minRad=None, raCol='RA', decCol='Dec')
Definition: matchFakes.py:146
lsst::daf::persistence.butler.Butler
Definition: butler.py:321
lsst.synpipe.matchFakes.getEllipse
def getEllipse(quad)
Definition: matchFakes.py:88
lsst::afw::geom::ellipses
Definition: Axes.h:42
lsst::afw::table._source.SourceCatalog
Definition: _source.py:33
lsst.synpipe.matchFakes.getMag
def getMag(flux, fluxerr, zeropoint)
Definition: matchFakes.py:81
lsst::afw::table::SchemaMapper
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:21
lsst::afw::table
Definition: table.dox:3
list
daf::base::PropertyList * list
Definition: fits.cc:913
type
table::Key< int > type
Definition: Detector.cc:163
lsst::pex::exceptions
Definition: Exception.h:37
lsst::daf::persistence
Definition: Utils.h:50
lsst::geom::Box2D
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
set
daf::base::PropertySet * set
Definition: fits.cc:912
lsst.synpipe.matchFakes.returnMatchTable
def returnMatchTable(rootDir, visit, ccdList, outfile=None, fakeCat=None, overwrite=False, filt=None, tol=1.0, pixMatch=False, multiband=False, reffMatch=False, pix=0.168, multijobs=1, includeMissing=True, minRad=None, raCol='RA', decCol='Dec')
Definition: matchFakes.py:511
lsst::afw::geom
Definition: frameSetUtils.h:40
Generated on Wed Jun 24 2020 18:10:12 for LSSTApplications by   doxygen 1.8.18