Coverage for python/lsst/meas/extensions/psfex/psfex.py: 0%

1# This file is part of meas_extensions_psfex.

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.

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/>.

22import os

23import re

25import numpy as np

26from astropy.io import fits

27try:

28 import matplotlib.pyplot as plt

29except ImportError:

30 plt = None

32import lsst.geom as geom

33import lsst.afw.geom as afwGeom

34from lsst.afw.fits import readMetadata

35import lsst.afw.image as afwImage

36import lsst.afw.table as afwTable

37import lsst.afw.display as afwDisplay

38from . import psfexLib

40afwDisplay.setDefaultMaskTransparency(75)

43def compute_fwhmrange(fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False)):

44 """Compute the FWHM range associated to a series of FWHM measurements.

45 AUTHOR E. Bertin (IAP, Leiden observatory & ESO)

46 VERSION 20/03/2008

48 Parameters

49 ----------

50 fwhm: iterable of `float`

51 Iterable of full width half-maxima.

52 maxvar: `float`

53 Maximum allowed FWHM variation.

54 minin: `float`

55 Minimum allowed FWHM.

56 maxin: `float`

57 Maximum allowed FWHM.

58 plot: `dict`, optional

59 Dict of plotting options.

61 Returns

62 -------

63 fmin: `float`

64 FWHM mode.

65 minout: `float`

66 Lower FWHM range.

67 maxout: `float`

68 Upper FWHM range.

69 """

71 nfwhm = len(fwhm)

72 fwhm.sort()

74 # Find the mode

75 nw = nfwhm//4

76 if nw < 4:

77 nw = 1

78 dfmin = psfexLib.BIG

79 fmin = 0.0

80 for i in range(nfwhm - nw):

81 df = fwhm[i + nw] - fwhm[i]

82 if df < dfmin:

83 dfmin = df

84 fmin = (fwhm[i + nw] + fwhm[i])/2.0

86 if nfwhm < 2:

87 fmin = fwhm[0]

89 dfmin = (maxvar + 1.0)**0.3333333

90 minout = fmin/dfmin if dfmin > 0.0 else 0.0

91 if minout < minin:

92 minout = minin

94 maxout = fmin*dfmin**2

95 if maxout > maxin:

96 maxout = maxin

98 if plt and plot.get("fwhmHistogram"):

99 plt.clf()

100 plt.hist(fwhm, nfwhm//10 + 1, normed=1, facecolor='g', alpha=0.75)

101 plt.xlabel("FWHM")

102 plt.axvline(fmin, color='red')

103 [plt.axvline(_, color='blue') for _ in (minout, maxout)]

104

105 input("Continue? ")

106

107 return fmin, minout, maxout

108

109

110def select_candidates(set, prefs, frmin, frmax,

111 flags, flux, fluxerr, rmsSize, elong, vignet,

112 plot=dict(), title=""):

113 maxbad = prefs.getBadpixNmax()

114 maxbadflag = prefs.getBadpixFlag()

115 maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100.0

116 minsn = prefs.getMinsn()

117

118 sn = flux/np.where(fluxerr > 0, fluxerr, 1)

119 sn[fluxerr <= 0] = -psfexLib.BIG

120 # ---- Apply some selection over flags, fluxes...

121 plotFlags = plot.get("showFlags") if plt else False

122 plotRejection = plot.get("showRejection") if plt else False

123

124 bad = flags & prefs.getFlagMask() != 0

125 set.setBadFlags(int(sum(bad != 0)))

126

127 if plotRejection:

128 selectionVectors = []

129 selectionVectors.append((bad, "flags %d" % sum(bad != 0)))

130

131 dbad = sn < minsn

132 set.setBadSN(int(sum(dbad)))

133 bad = np.logical_or(bad, dbad)

134 if plotRejection:

135 selectionVectors.append((dbad, "S/N %d" % sum(dbad)))

136

137 dbad = rmsSize < frmin

138 set.setBadFrmin(int(sum(dbad)))

139 bad = np.logical_or(bad, dbad)

140 if plotRejection:

141 selectionVectors.append((dbad, "frmin %d" % sum(dbad)))

142

143 dbad = rmsSize > frmax

144 set.setBadFrmax(int(sum(dbad)))

145 bad = np.logical_or(bad, dbad)

146 if plotRejection:

147 selectionVectors.append((dbad, "frmax %d" % sum(dbad)))

148

149 dbad = elong > maxelong

150 set.setBadElong(int(sum(dbad)))

151 bad = np.logical_or(bad, dbad)

152 if plotRejection:

153 selectionVectors.append((dbad, "elong %d" % sum(dbad)))

154

155 # -- ... and check the integrity of the sample

156 if maxbadflag:

157 nbad = np.array([(v <= -psfexLib.BIG).sum() for v in vignet])

158 dbad = nbad > maxbad

159 set.setBadPix(int(sum(dbad)))

160 bad = np.logical_or(bad, dbad)

161 if plotRejection:

162 selectionVectors.append((dbad, "badpix %d" % sum(dbad)))

163

164 good = np.logical_not(bad)

165 if plotFlags or plotRejection:

166 imag = -2.5*np.log10(flux)

167 plt.clf()

168

169 alpha = 0.5

170 if plotFlags:

171 labels = getFlags()

172

173 isSet = np.where(flags == 0x0)[0]

174 plt.plot(imag[isSet], rmsSize[isSet], 'o', alpha=alpha, label="good")

175

176 for i in range(16):

177 mask = 1 << i

178 if mask & prefs.getFlagMask():

179 isSet = np.where(np.bitwise_and(flags, mask))[0]

180 if isSet.any():

181 plt.plot(imag[isSet], rmsSize[isSet], 'o', alpha=alpha, label=labels[mask])

182 else:

183 for bad, label in selectionVectors:

184 plt.plot(imag[bad], rmsSize[bad], 'o', alpha=alpha, label=label)

185

186 plt.plot(imag[good], rmsSize[good], 'o', color="black", label="selected")

187 [plt.axhline(_, color='red') for _ in [frmin, frmax]]

188 plt.xlim(np.median(imag[good]) + 5*np.array([-1, 1]))

189 plt.ylim(-0.1, 2*frmax)

190 plt.legend(loc=2)

191 plt.xlabel("Instrumental Magnitude")

192 plt.ylabel("rmsSize")

193 plt.title("%s %d selected" % (title, sum(good)))

194

195 input("Continue? ")

196

197 return good

198

199

200def showPsf(psf, set, ext=None, wcsData=None, trim=0, nspot=5,

201 diagnostics=False, outDir="", frame=None, title=None):

202 """Show a PSF on display (e.g., ds9)

203 """

204

205 if ext is not None:

206 psf = psf[ext]

207

208 if wcsData:

209 if ext is not None:

210 wcsData = wcsData[ext]

211 wcs, naxis1, naxis2 = wcsData

212 else:

213 wcs, naxis1, naxis2 = None, None, None

214

215 naxis = [naxis1, naxis2]

216 for i in range(2):

217 if naxis[i] is None:

218 # cmin, cmax are the range of input star positions

219 cmin, cmax = [set.getContextOffset(i) + d*set.getContextScale(i) for d in (-0.5, 0.5)]

220 naxis[i] = cmax + cmin # a decent guess

221

222 if naxis[0] > naxis[1]:

223 nx, ny = int(nspot*naxis[0]/float(naxis[1]) + 0.5), nspot

224 else:

225 nx, ny = nspot, int(nspot*naxis[1]/float(naxis[0]) + 0.5)

226

227 mos = afwDisplay.utils.Mosaic(gutter=2, background=0.02)

228

229 xpos, ypos = np.linspace(0, naxis[0], nx), np.linspace(0, naxis[1], ny)

230 for y in ypos:

231 for x in xpos:

232 psf.build(x, y)

233

234 im = afwImage.ImageF(*psf.getLoc().shape)

235 im.getArray()[:] = psf.getLoc()

236 im /= float(im.getArray().max())

237 if trim:

238 if trim > im.getHeight()//2:

239 trim = im.getHeight()//2

240

241 im = im[trim:-trim, trim:-trim]

242

243 mos.append(im)

244

245 mosaic = mos.makeMosaic(mode=nx)

246 if frame is not None:

247 afwDisplay.Display(frame=frame).mtv(mosaic, title=title)

248 #

249 # Figure out the WCS for the mosaic

250 #

251 pos = []

252 pos.append([geom.PointD(0, 0), wcs.pixelToSky(geom.PointD(0, 0))])

253 pos.append([geom.PointD(*mosaic.getDimensions()), wcs.pixelToSky(geom.PointD(naxis1, naxis2))])

254

255 CD = []

256 for i in range(2):

257 delta = pos[1][1][i].asDegrees() - pos[0][1][i].asDegrees()

258 CD.append(delta/(pos[1][0][i] - pos[0][0][i]))

259 CD = np.array(CD)

260 CD.shape = (2, 2)

261 mosWcs = afwGeom.makeSkyWcs(crval=pos[0][0], crpix=pos[0][1], cdMatrix=CD)

262

263 if ext is not None:

264 title = "%s-%d" % (title, ext)

265

266 if frame is not None:

267 afwDisplay.Display(frame=frame).mtv(mosaic, title=title, wcs=mosWcs)

268

269 if diagnostics:

270 outFile = "%s-mod.fits" % title

271 if outDir:

272 outFile = os.path.join(outDir, outFile)

273 mosaic.writeFits(outFile, mosWcs.getFitsMetadata())

274

275 mos = afwDisplay.utils.Mosaic(gutter=4, background=0.002)

276 for i in range(set.getNsample()):

277 s = set.getSample(i)

278 if ext is not None and s.getExtindex() != ext:

279 continue

280

281 smos = afwDisplay.utils.Mosaic(gutter=2, background=-0.003)

282 for func in [s.getVig, s.getVigResi]:

283 arr = func()

284 if func == s.getVig:

285 norm = float(arr.max()) if True else s.getNorm()

286

287 arr /= norm

288 im = afwImage.ImageF(*arr.shape)

289 im.getArray()[:] = arr

290 smos.append(im)

291

292 mos.append(smos.makeMosaic(mode="x"))

293

294 mosaic = mos.makeMosaic(title=title)

295

296 if frame is not None:

297 afwDisplay.Display(frame=frame + 1).mtv(mosaic, title=title)

298

299 if diagnostics:

300 outFile = "%s-psfstars.fits" % title

301 if outDir:

302 outFile = os.path.join(outDir, outFile)

303

304 mosaic.writeFits(outFile)

305

306

307def getFlags(tab=None):

308 flagKeys = [

309 "base_PixelFlags_flag_edge",

310 "base_PixelFlags_flag_nodata",

311 # "base_PixelFlags_flag_interpolated",

312 # "base_PixelFlags_flag_interpolatedCenter",

313 # "base_PixelFlags_flag_saturated",

314 "base_PixelFlags_flag_saturatedCenter",

315 # "base_PixelFlags_flag_cr",

316 "base_PixelFlags_flag_crCenter",

317 "base_PixelFlags_flag_bad",

318 "base_PsfFlux_flag",

319 "parent",

320 ]

321

322 if tab is None:

323 flags = {}

324 for i, k in enumerate(flagKeys):

325 flags[1 << i] = re.sub(r"\_flag", "",

326 re.sub(r"^base\_", "", re.sub(r"^base\_PixelFlags\_flag\_", "", k)))

327 else:

328 flags = 0

329 for i, k in enumerate(flagKeys):

330 if k == "parent":

331 try:

332 isSet = tab.get("deblend_nChild") > 0

333 except KeyError:

334 isSet = 0

335 else:

336 isSet = tab.get(k)

337 flags = np.bitwise_or(flags, np.where(isSet, 1 << i, 0))

338

339 return flags

340

341

342def guessCalexp(fileName):

343 for guess in [

344 re.sub("/src", r"", fileName),

345 re.sub("(SRC([^.]+))", r"CORR\2-exp", fileName),

346 ]:

347 if guess != fileName and os.path.exists(guess):

348 return guess

349

350 raise RuntimeError("Unable to find a calexp to go with %s" % fileName)

351

352

353def makeitLsst(prefs, context, saveWcs=False, plot=dict()):

354 """This is the python wrapper that reads lsst tables

355 """

356 # Create an array of PSFs (one PSF for each extension)

357 if prefs.getVerboseType() != prefs.QUIET:

358 print("----- %d input catalogues:" % prefs.getNcat())

359

360 if saveWcs: # only needed for making plots

361 wcssList = []

362

363 fields = psfexLib.vectorField()

364 for cat in prefs.getCatalogs():

365 field = psfexLib.Field(cat)

366 wcss = []

367 wcssList.append(wcss)

368 with fits.open(cat):

369 # Hack: I want the WCS so I'll guess where the calexp is to be

370 # found

371 calexpFile = guessCalexp(cat)

372 md = readMetadata(calexpFile)

373 wcs = afwGeom.makeSkyWcs(md)

374

375 if not wcs:

376 cdMatrix = np.array([1.0, 0.0, 0.0, 1.0])

377 cdMatrix.shape = (2, 2)

378 wcs = afwGeom.makeSkyWcs(crpix=geom.PointD(0, 0),

379 crval=geom.SpherePoint(0.0, 0.0, geom.degrees),

380 cdMatrix=cdMatrix)

381

382 naxis1, naxis2 = md.getScalar("NAXIS1"), md.getScalar("NAXIS2")

383 # Find how many rows there are in the catalogue

384 md = readMetadata(cat)

385

386 field.addExt(wcs, naxis1, naxis2, md.getScalar("NAXIS2"))

387 if saveWcs:

388 wcss.append((wcs, naxis1, naxis2))

389

390 field.finalize()

391 fields.append(field)

392

393 fields[0].getNext() # number of extensions

394

395 prefs.getPsfStep()

396

397 sets = psfexLib.vectorSet()

398 for set in load_samplesLsst(prefs, context, plot=plot):

399 sets.append(set)

400

401 psfexLib.makeit(fields, sets)

402

403 ret = [[f.getPsfs() for f in fields], sets]

404 if saveWcs:

405 ret.append(wcssList)

406

407 return ret

408

409

410def read_samplesLsst(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, nobj,

411 plot=dict(showFlags=False, showRejection=False)):

412 # allocate a new set iff set is None

413 if not set:

414 set = psfexLib.Set(context)

415

416 cmin, cmax = None, None

417 if set.getNcontext():

418 cmin = np.empty(set.getNcontext())

419 cmax = np.empty(set.getNcontext())

420 for i in range(set.getNcontext()):

421 if set.getNsample():

422 cmin[i] = set.getContextOffset(i) - set.getContextScale(i)/2.0

423 cmax[i] = cmin[i] + set.getContextScale(i)

424 else:

425 cmin[i] = psfexLib.BIG

426 cmax[i] = -psfexLib.BIG

427 #

428 # Read data

429 #

430 tab = afwTable.SourceCatalog.readFits(filename)

431

432 centroid = tab.getCentroidDefinition()

433 xm = tab.get("%s.x" % centroid)

434 ym = tab.get("%s.y" % centroid)

435

436 shape = tab.getShapeDefinition()

437 ixx = tab.get("%s.xx" % shape)

438 iyy = tab.get("%s.yy" % shape)

439

440 rmsSize = np.sqrt(0.5*(ixx + iyy))

441 elong = 0.5*(ixx - iyy)/(ixx + iyy)

442

443 flux = tab.get(prefs.getPhotfluxRkey())

444 fluxErr = tab.get(prefs.getPhotfluxerrRkey())

445 flags = getFlags(tab)

446

447 #

448 # Now the VIGNET data

449 #

450 vigw, vigh = 35, 35 # [prefs.getPsfsize()[i] for i in range(2)]

451 if set.empty():

452 set.setVigSize(vigw, vigh)

453

454 vignet = np.empty(nobj*vigw*vigh, "float32").reshape(nobj, vigw, vigh)

455

456 # Hack: I want the WCS so I'll guess where the calexp is to be found

457 calexpFile = guessCalexp(filename)

458 mi = afwImage.MaskedImageF(calexpFile)

459 backnoise2 = np.median(mi.getVariance().getArray())

460 gain = 1.0

461

462 edgeBit = [k for k, v in getFlags().items() if v == "edge"][0]

463

464 for i, xc, yc in zip(range(nobj), xm, ym):

465 try:

466 x, y = int(xc), int(yc)

467 except ValueError:

468 flags[i] |= edgeBit # mark star as bad

469

470 try:

471 pstamp = mi[x - vigw//2:x + vigw//2 + 1, y - vigh//2:y + vigh//2 + 1]

472 vignet[i] = pstamp.getImage().getArray().transpose()

473 except Exception:

474 flags[i] |= edgeBit # mark star as bad

475

476 # Try to load the set of context keys

477 pc = 0

478 contextvalp = []

479 for i, key in enumerate(context.getName()):

480 if context.getPcflag(i):

481 contextvalp.append(pcval[pc])

482 pc += 1

483 elif key[0] == ':':

484 try:

485 contextvalp.append(tab.header[key[1:]])

486 except KeyError:

487 raise RuntimeError("*Error*: %s parameter not found in the header of %s" %

488 (key[1:], filename))

489 else:

490 try:

491 contextvalp.append(tab.get(key))

492 except KeyError:

493 raise RuntimeError("*Error*: %s parameter not found in the header of %s" %

494 (key, filename))

495 set.setContextname(i, key)

496

497 # Now examine each vector of the shipment

498 good = select_candidates(set, prefs, frmin, frmax,

499 flags, flux, fluxErr, rmsSize, elong, vignet,

500 plot=plot, title="%s[%d]" % (filename, ext + 1) if next > 1 else filename)

501 #

502 # Insert the good candidates into the set

503 #

504 if not vignet.dtype.isnative:

505 # without the swap setVig fails with

506 # "ValueError: 'unaligned arrays cannot be converted to C++'"

507 vignet = vignet.byteswap()

508

509 for i in np.where(good)[0]:

510 sample = set.newSample()

511 sample.setCatindex(catindex)

512 sample.setExtindex(ext)

513

514 sample.setVig(vignet[i])

515 sample.setNorm(float(flux[i]))

516 sample.setBacknoise2(backnoise2)

517 sample.setGain(gain)

518 sample.setX(float(xm[i]))

519 sample.setY(float(ym[i]))

520 sample.setFluxrad(float(rmsSize[i]))

521

522 for j in range(set.getNcontext()):

523 sample.setContext(j, float(contextvalp[j][i]))

524

525 set.finiSample(sample, prefs.getProfAccuracy())

526

527 # ---- Update min and max

528 for j in range(set.getNcontext()):

529 cmin[j] = contextvalp[j][good].min()

530 cmax[j] = contextvalp[j][good].max()

531

532 # Update the scaling

533 if set.getNsample():

534 for i in range(set.getNcontext()):

535 set.setContextScale(i, cmax[i] - cmin[i])

536 set.setContextOffset(i, (cmin[i] + cmax[i])/2.0)

537

538 # Don't waste memory!

539 set.trimMemory()

540

541 return set

542

543

544def load_samplesLsst(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict()):

545 minsn = prefs.getMinsn()

546 maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100

547 min = prefs.getFwhmrange()[0]

548 max = prefs.getFwhmrange()[1]

549

550 filenames = prefs.getCatalogs()

551

552 ncat = len(filenames)

553 fwhmmin = np.empty(ncat)

554 fwhmmax = np.empty(ncat)

555

556 if not prefs.getAutoselectFlag():

557 fwhmmin = np.zeros(ncat) + prefs.getFwhmrange()[0]

558 fwhmmax = np.zeros(ncat) + prefs.getFwhmrange()[1]

559 fwhmmode = (fwhmmin + fwhmmax)/2.0

560 else:

561 fwhms = {}

562

563 # -- Try to estimate the most appropriate Half-light Radius range

564 # -- Get the Half-light radii

565 nobj = 0

566 for i, fileName in enumerate(filenames):

567 fwhms[i] = []

568

569 if prefs.getVerboseType() != prefs.QUIET:

570 print("Examining Catalog #%d" % (i+1))

571

572 # ---- Read input catalog

573 tab = afwTable.SourceCatalog.readFits(fileName)

574

575 # -------- Fill the FWHM array

576 shape = tab.getShapeDefinition()

577 ixx = tab.get("%s.xx" % shape)

578 iyy = tab.get("%s.yy" % shape)

579

580 rmsSize = np.sqrt(0.5*(ixx + iyy))

581 elong = 0.5*(ixx - iyy)/(ixx + iyy)

582

583 flux = tab.get(prefs.getPhotfluxRkey())

584 fluxErr = tab.get(prefs.getPhotfluxerrRkey())

585

586 flags = getFlags(tab)

587

588 good = np.logical_and(flux/fluxErr > minsn,

589 np.logical_not(flags & prefs.getFlagMask()))

590 good = np.logical_and(good, elong < maxelong)

591 fwhm = 2.0*rmsSize

592 good = np.logical_and(good, fwhm >= min)

593 good = np.logical_and(good, fwhm < max)

594 fwhms[i] = fwhm[good]

595

596 if prefs.getVarType() == prefs.VAR_NONE:

597 if nobj:

598 fwhms_all = np.empty(sum([len(f) for f in fwhms.values()]))

599 i = 0

600 for f in fwhms.values():

601 fwhms_all[i:len(f)] = f

602 i += len(f)

603 mode, min, max = compute_fwhmrange(fwhms_all, prefs.getMaxvar(),

604 prefs.getFwhmrange()[0], prefs.getFwhmrange()[1],

605 plot=plot)

606 else:

607 raise RuntimeError("No source with appropriate FWHM found!!")

608 mode = min = max = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)

609

610 fwhmmin = np.zeros(ncat) + min

611 fwhmmax = np.zeros(ncat) + max

612 fwhmmode = np.zeros(ncat) + mode

613 else:

614 fwhmmode = np.empty(ncat)

615 fwhmmin = np.empty(ncat)

616 fwhmmax = np.empty(ncat)

617

618 for i in range(ncat):

619 nobj = len(fwhms[i])

620 if (nobj):

621 fwhmmode[i], fwhmmin[i], fwhmmax[i] = \

622 compute_fwhmrange(fwhms[i], prefs.getMaxvar(),

623 prefs.getFwhmrange()[0], prefs.getFwhmrange()[1], plot=plot)

624 else:

625 raise RuntimeError("No source with appropriate FWHM found!!")

626 fwhmmode[i] = fwhmmin[i] = fwhmmax[i] = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)

627

628 # Read the samples

629 mode = psfexLib.BIG # mode of FWHM distribution

630

631 sets = []

632 for i, fileName in enumerate(filenames):

633 set = None

634 for ext in range(next):

635 set = read_samplesLsst(prefs, set, fileName, fwhmmin[i]/2.0, fwhmmax[i]/2.0,

636 ext, next, i, context,

637 context.getPc(i) if context.getNpc() else None, nobj, plot=plot)

638

639 if fwhmmode[i] < mode:

640 mode = fwhmmode[i]

641

642 set.setFwhm(mode)

643

644 if prefs.getVerboseType() != prefs.QUIET:

645 if set.getNsample():

646 print("%d samples loaded." % set.getNsample())

647 else:

648 raise RuntimeError("No appropriate source found!!")

649

650 sets.append(set)

651

652 return sets

Coverage for python / lsst / meas / extensions / psfex / psfex.py: 0%

396 statements