219 raCol=
'RA', decCol=
'Dec'):
221 Get list of sources which agree in pixel position with fake ones with tol. 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 228 The outputs can include extraCols as long as they are one of: 229 zeropoint, visit, ccd, thetaNorth, pixelScale 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 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 239 coaddData =
"deepCoadd_calexp" 240 coaddMeta =
"deepCoadd_calexp_md" 242 availExtras = {
'zeropoint': {
'type': float,
'doc':
'zeropoint'},
243 'visit': {
'type': int,
'doc':
'visit id'},
244 'ccd': {
'type': int,
'doc':
'ccd id'},
246 'doc':
'angle to north'},
247 'pixelScale': {
'type': float,
248 'doc':
'pixelscale in arcsec/pixel'}}
250 if not np.in1d(extraCols,
list(availExtras.keys())).
all():
251 print(
"extraCols must be in ", availExtras)
254 if 'filter' not in dataId:
255 sources = butler.get(
'src', dataId,
256 flags=lsst.afw.table.SOURCE_IO_NO_FOOTPRINTS,
258 cal = butler.get(
'calexp', dataId, immediate=
True)
259 cal_md = butler.get(
'calexp_md', dataId, immediate=
True)
261 meas = butler.get(
'deepCoadd_meas', dataId,
264 force = butler.get(
'deepCoadd_forced_src', dataId,
268 cal = butler.get(coaddData, dataId, immediate=
True)
269 cal_md = butler.get(coaddMeta, dataId, immediate=
True)
271 print(
"skipping", dataId)
274 if (
'pixelScale' in extraCols)
or (
'thetaNorth' in extraCols):
276 availExtras[
'pixelScale'][
'value'] = wcs.getPixelScale().asArcseconds()
280 xMid = cal.getWidth() // 2
281 yMid = cal.getHeight() // 2
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
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
297 if radecMatch
is None:
300 if minRad
is not None:
301 print(
"# The min matching radius is %4.1f pixel" % minRad)
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?')
336 for extraName
in set(extraCols).intersection(availExtras):
337 newSchema.addField(extraName, type=availExtras[extraName][
'type'],
338 doc=availExtras[extraName][
'doc'])
341 srcList.reserve(sum([len(s)
for s
in srcIndex.values()]) +
342 (0
if not includeMissing
else list(srcIndex.values()).count([])))
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:
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 356 if includeMissing
and (nMatched == 0):
357 newRec = srcList.addNew()
358 newRec.set(
'fakeId', ident)
360 newRec.set(
'nMatched', 0)
361 newRec.set(
'rMatched', rMatched)
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() -
372 newRec.set(
'fakeOffX', offsetX)
373 offsetY = (sources[
int(ss)].get(centroidKey).getY() -
375 newRec.set(
'fakeOffY', offsetY)
376 newRec.set(
'fakeOffR', np.sqrt(offsetX ** 2.0 + offsetY ** 2.0))
378 if int(ss) ==
int(srcClose[ident]):
379 newRec.set(
'fakeClosest',
True)
381 newRec.set(
'fakeClosest',
False)
384 srcList = srcList.copy(deep=
True)
386 for extraName
in set(extraCols).intersection(availExtras):
387 tempCol = srcList.get(extraName)
388 tempCol.fill(availExtras[extraName][
'value'])
A floating-point coordinate rectangle geometry.
daf::base::PropertySet * set
A class representing an angle.
bool all(CoordinateExpr< N > const &expr) noexcept
Return true if all elements are true.
def combineWithForce(meas, force)
def getFakeMatchesHeader(cal_md, sources, tol=1.0)
SortedCatalogT< SourceRecord > SourceCatalog
def getFakeMatchesRaDec(sources, radecCatFile, bbox, wcs, tol=1.0, reffMatch=False, pix=0.168, minRad=None, raCol='RA', decCol='Dec')
daf::base::PropertyList * list