Coverage for python/lsst/pipe/tasks/functors.py: 39%

1# This file is part of pipe_tasks.

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

22__all__ = ["init_fromDict", "Functor", "CompositeFunctor", "mag_aware_eval",

23 "CustomFunctor", "Column", "Index", "CoordColumn", "RAColumn",

24 "DecColumn", "SinglePrecisionFloatColumn", "HtmIndex20", "fluxName", "fluxErrName", "Mag",

25 "MagErr", "MagDiff", "Color", "DeconvolvedMoments", "SdssTraceSize",

26 "PsfSdssTraceSizeDiff", "HsmTraceSize", "PsfHsmTraceSizeDiff",

27 "HsmFwhm", "E1", "E2", "RadiusFromQuadrupole", "LocalWcs",

28 "ComputePixelScale", "ConvertPixelToArcseconds",

29 "ConvertPixelSqToArcsecondsSq",

30 "ConvertDetectorAngleToPositionAngle",

31 "ConvertDetectorAngleErrToPositionAngleErr",

32 "ReferenceBand", "Photometry",

33 "NanoJansky", "NanoJanskyErr", "LocalPhotometry", "LocalNanojansky",

34 "LocalNanojanskyErr", "LocalDipoleMeanFlux",

35 "LocalDipoleMeanFluxErr", "LocalDipoleDiffFlux",

36 "LocalDipoleDiffFluxErr", "Ebv",

37 "MomentsIuuSky", "MomentsIvvSky", "MomentsIuvSky",

38 "CorrelationIuuSky", "CorrelationIvvSky", "CorrelationIuvSky",

39 "PositionAngleFromMoments", "PositionAngleFromCorrelation",

40 "SemimajorAxisFromMoments", "SemimajorAxisFromCorrelation",

41 "SemiminorAxisFromMoments", "SemiminorAxisFromCorrelation",

42 ]

44import logging

45import os

46import os.path

47import re

48import warnings

49from contextlib import redirect_stdout

50from itertools import product

52import astropy.units as u

53import lsst.geom as geom

54import lsst.sphgeom as sphgeom

55import numpy as np

56import pandas as pd

57import yaml

58from astropy.coordinates import SkyCoord

59from lsst.daf.butler import DeferredDatasetHandle

60from lsst.pipe.base import InMemoryDatasetHandle

61from lsst.utils import doImport

62from lsst.utils.introspection import get_full_type_name

65def init_fromDict(initDict, basePath='lsst.pipe.tasks.functors',

66 typeKey='functor', name=None):

67 """Initialize an object defined in a dictionary.

69 The object needs to be importable as f'{basePath}.{initDict[typeKey]}'.

70 The positional and keyword arguments (if any) are contained in "args" and

71 "kwargs" entries in the dictionary, respectively.

72 This is used in `~lsst.pipe.tasks.functors.CompositeFunctor.from_yaml` to

73 initialize a composite functor from a specification in a YAML file.

75 Parameters

76 ----------

77 initDict : dictionary

78 Dictionary describing object's initialization.

79 Must contain an entry keyed by ``typeKey`` that is the name of the

80 object, relative to ``basePath``.

81 basePath : str

82 Path relative to module in which ``initDict[typeKey]`` is defined.

83 typeKey : str

84 Key of ``initDict`` that is the name of the object (relative to

85 ``basePath``).

86 """

87 initDict = initDict.copy()

88 # TO DO: DM-21956 We should be able to define functors outside this module

89 pythonType = doImport(f'{basePath}.{initDict.pop(typeKey)}')

90 args = []

91 if 'args' in initDict:

92 args = initDict.pop('args')

93 if isinstance(args, str):

94 args = [args]

95 try:

96 element = pythonType(*args, **initDict)

97 except Exception as e:

98 message = f'Error in constructing functor "{name}" of type {pythonType.__name__} with args: {args}'

99 raise type(e)(message, e.args)

100 return element

101

102

103class Functor(object):

104 """Define and execute a calculation on a DataFrame or Handle holding a

105 DataFrame.

106

107 The `__call__` method accepts either a `~pandas.DataFrame` object or a

108 `~lsst.daf.butler.DeferredDatasetHandle` or

109 `~lsst.pipe.base.InMemoryDatasetHandle`, and returns the

110 result of the calculation as a single column.

111 Each functor defines what columns are needed for the calculation, and only

112 these columns are read from the dataset handle.

113

114 The action of `__call__` consists of two steps: first, loading the

115 necessary columns from disk into memory as a `~pandas.DataFrame` object;

116 and second, performing the computation on this DataFrame and returning the

117 result.

118

119 To define a new `Functor`, a subclass must define a `_func` method,

120 that takes a `~pandas.DataFrame` and returns result in a `~pandas.Series`.

121 In addition, it must define the following attributes:

122

123 * `_columns`: The columns necessary to perform the calculation

124 * `name`: A name appropriate for a figure axis label

125 * `shortname`: A name appropriate for use as a dictionary key

126

127 On initialization, a `Functor` should declare what band (``filt`` kwarg)

128 and dataset (e.g. ``'ref'``, ``'meas'``, ``'forced_src'``) it is intended

129 to be applied to.

130 This enables the `_get_data` method to extract the proper columns from the

131 underlying data.

132 If not specified, the dataset will fall back on the `_defaultDataset`

133 attribute.

134 If band is not specified and ``dataset`` is anything other than ``'ref'``,

135 then an error will be raised when trying to perform the calculation.

136

137 Originally, `Functor` was set up to expect datasets formatted like the

138 ``deepCoadd_obj`` dataset; that is, a DataFrame with a multi-level column

139 index, with the levels of the column index being ``band``, ``dataset``, and

140 ``column``.

141 It has since been generalized to apply to DataFrames without multi-level

142 indices and multi-level indices with just ``dataset`` and ``column``

143 levels.

144 In addition, the `_get_data` method that reads the columns from the

145 underlying data will return a DataFrame with column index levels defined by

146 the `_dfLevels` attribute; by default, this is ``column``.

147

148 The `_dfLevels` attributes should generally not need to be changed, unless

149 `_func` needs columns from multiple filters or datasets to do the

150 calculation.

151 An example of this is the `~lsst.pipe.tasks.functors.Color` functor, for

152 which `_dfLevels = ('band', 'column')`, and `_func` expects the DataFrame

153 it gets to have those levels in the column index.

154

155 Parameters

156 ----------

157 filt : str

158 Band upon which to do the calculation.

159

160 dataset : str

161 Dataset upon which to do the calculation (e.g., 'ref', 'meas',

162 'forced_src').

163 """

164

165 _defaultDataset = 'ref'

166 _dfLevels = ('column',)

167 _defaultNoDup = False

168

169 def __init__(self, filt=None, dataset=None, noDup=None):

170 self.filt = filt

171 self.dataset = dataset if dataset is not None else self._defaultDataset

172 self._noDup = noDup

173 self.log = logging.getLogger(type(self).__name__)

174

175 @property

176 def noDup(self):

177 """Do not explode by band if used on object table."""

178 if self._noDup is not None:

179 return self._noDup

180 else:

181 return self._defaultNoDup

182

183 @property

184 def columns(self):

185 """Columns required to perform calculation."""

186 if not hasattr(self, '_columns'):

187 raise NotImplementedError('Must define columns property or _columns attribute')

188 return self._columns

189

190 def _get_data_columnLevels(self, data, columnIndex=None):

191 """Gets the names of the column index levels.

192

193 This should only be called in the context of a multilevel table.

194

195 Parameters

196 ----------

197 data : various

198 The data to be read, can be a

199 `~lsst.daf.butler.DeferredDatasetHandle` or

200 `~lsst.pipe.base.InMemoryDatasetHandle`.

201 columnIndex (optional): pandas `~pandas.Index` object

202 If not passed, then it is read from the

203 `~lsst.daf.butler.DeferredDatasetHandle`

204 for `~lsst.pipe.base.InMemoryDatasetHandle`.

205 """

206 if columnIndex is None:

207 columnIndex = data.get(component="columns")

208 return columnIndex.names

209

210 def _get_data_columnLevelNames(self, data, columnIndex=None):

211 """Gets the content of each of the column levels for a multilevel

212 table.

213 """

214 if columnIndex is None:

215 columnIndex = data.get(component="columns")

216

217 columnLevels = columnIndex.names

218 columnLevelNames = {

219 level: list(np.unique(np.array([c for c in columnIndex])[:, i]))

220 for i, level in enumerate(columnLevels)

221 }

222 return columnLevelNames

223

224 def _colsFromDict(self, colDict, columnIndex=None):

225 """Converts dictionary column specficiation to a list of columns."""

226 new_colDict = {}

227 columnLevels = self._get_data_columnLevels(None, columnIndex=columnIndex)

228

229 for i, lev in enumerate(columnLevels):

230 if lev in colDict:

231 if isinstance(colDict[lev], str):

232 new_colDict[lev] = [colDict[lev]]

233 else:

234 new_colDict[lev] = colDict[lev]

235 else:

236 new_colDict[lev] = columnIndex.levels[i]

237

238 levelCols = [new_colDict[lev] for lev in columnLevels]

239 cols = list(product(*levelCols))

240 colsAvailable = [col for col in cols if col in columnIndex]

241 return colsAvailable

242

243 def multilevelColumns(self, data, columnIndex=None, returnTuple=False):

244 """Returns columns needed by functor from multilevel dataset.

245

246 To access tables with multilevel column structure, the

247 `~lsst.daf.butler.DeferredDatasetHandle` or

248 `~lsst.pipe.base.InMemoryDatasetHandle` needs to be passed

249 either a list of tuples or a dictionary.

250

251 Parameters

252 ----------

253 data : various

254 The data as either `~lsst.daf.butler.DeferredDatasetHandle`, or

255 `~lsst.pipe.base.InMemoryDatasetHandle`.

256 columnIndex (optional): pandas `~pandas.Index` object

257 Either passed or read in from

258 `~lsst.daf.butler.DeferredDatasetHandle`.

259 `returnTuple` : `bool`

260 If true, then return a list of tuples rather than the column

261 dictionary specification.

262 This is set to `True` by `CompositeFunctor` in order to be able to

263 combine columns from the various component functors.

264

265 """

266 if not isinstance(data, (DeferredDatasetHandle, InMemoryDatasetHandle)):

267 raise RuntimeError(f"Unexpected data type. Got {get_full_type_name(data)}.")

268

269 if columnIndex is None:

270 columnIndex = data.get(component="columns")

271

272 # Confirm that the dataset has the column levels the functor is

273 # expecting it to have.

274 columnLevels = self._get_data_columnLevels(data, columnIndex)

275

276 columnDict = {'column': self.columns,

277 'dataset': self.dataset}

278 if self.filt is None:

279 columnLevelNames = self._get_data_columnLevelNames(data, columnIndex)

280 if "band" in columnLevels:

281 if self.dataset == "ref":

282 columnDict["band"] = columnLevelNames["band"][0]

283 else:

284 raise ValueError(f"'filt' not set for functor {self.name}"

285 f"(dataset {self.dataset}) "

286 "and DataFrame "

287 "contains multiple filters in column index. "

288 "Set 'filt' or set 'dataset' to 'ref'.")

289 else:

290 columnDict['band'] = self.filt

291

292 if returnTuple:

293 return self._colsFromDict(columnDict, columnIndex=columnIndex)

294 else:

295 return columnDict

296

297 def _func(self, df, dropna=True):

298 raise NotImplementedError('Must define calculation on DataFrame')

299

300 def _get_columnIndex(self, data):

301 """Return columnIndex."""

302

303 if isinstance(data, (DeferredDatasetHandle, InMemoryDatasetHandle)):

304 return data.get(component="columns")

305 else:

306 return None

307

308 def _get_data(self, data):

309 """Retrieve DataFrame necessary for calculation.

310

311 The data argument can be a `~pandas.DataFrame`, a

312 `~lsst.daf.butler.DeferredDatasetHandle`, or

313 an `~lsst.pipe.base.InMemoryDatasetHandle`.

314

315 Returns a DataFrame upon which `self._func` can act.

316 """

317 # We wrap a DataFrame in a handle here to take advantage of the

318 # DataFrame delegate DataFrame column wrangling abilities.

319 if isinstance(data, pd.DataFrame):

320 _data = InMemoryDatasetHandle(data, storageClass="DataFrame")

321 elif isinstance(data, (DeferredDatasetHandle, InMemoryDatasetHandle)):

322 _data = data

323 else:

324 raise RuntimeError(f"Unexpected type provided for data. Got {get_full_type_name(data)}.")

325

326 # First thing to do: check to see if the data source has a multilevel

327 # column index or not.

328 columnIndex = self._get_columnIndex(_data)

329 is_multiLevel = isinstance(columnIndex, pd.MultiIndex)

330

331 # Get proper columns specification for this functor.

332 if is_multiLevel:

333 columns = self.multilevelColumns(_data, columnIndex=columnIndex)

334 else:

335 columns = self.columns

336

337 # Load in-memory DataFrame with appropriate columns the gen3 way.

338 df = _data.get(parameters={"columns": columns})

339

340 # Drop unnecessary column levels.

341 if is_multiLevel:

342 df = self._setLevels(df)

343

344 return df

345

346 def _setLevels(self, df):

347 levelsToDrop = [n for n in df.columns.names if n not in self._dfLevels]

348 df.columns = df.columns.droplevel(levelsToDrop)

349 return df

350

351 def _dropna(self, vals):

352 return vals.dropna()

353

354 def __call__(self, data, dropna=False):

355 df = self._get_data(data)

356 try:

357 vals = self._func(df)

358 except Exception as e:

359 self.log.error("Exception in %s call: %s: %s", self.name, type(e).__name__, e)

360 vals = self.fail(df)

361 if dropna:

362 vals = self._dropna(vals)

363

364 return vals

365

366 def difference(self, data1, data2, **kwargs):

367 """Computes difference between functor called on two different

368 DataFrame/Handle objects.

369 """

370 return self(data1, **kwargs) - self(data2, **kwargs)

371

372 def fail(self, df):

373 return pd.Series(np.full(len(df), np.nan), index=df.index)

374

375 @property

376 def name(self):

377 """Full name of functor (suitable for figure labels)."""

378 return NotImplementedError

379

380 @property

381 def shortname(self):

382 """Short name of functor (suitable for column name/dict key)."""

383 return self.name

384

385

386class CompositeFunctor(Functor):

387 """Perform multiple calculations at once on a catalog.

388

389 The role of a `CompositeFunctor` is to group together computations from

390 multiple functors.

391 Instead of returning `~pandas.Series` a `CompositeFunctor` returns a

392 `~pandas.DataFrame`, with the column names being the keys of ``funcDict``.

393

394 The `columns` attribute of a `CompositeFunctor` is the union of all columns

395 in all the component functors.

396

397 A `CompositeFunctor` does not use a `_func` method itself; rather, when a

398 `CompositeFunctor` is called, all its columns are loaded at once, and the

399 resulting DataFrame is passed to the `_func` method of each component

400 functor.

401 This has the advantage of only doing I/O (reading from parquet file) once,

402 and works because each individual `_func` method of each component functor

403 does not care if there are *extra* columns in the DataFrame being passed;

404 only that it must contain *at least* the `columns` it expects.

405

406 An important and useful class method is `from_yaml`, which takes as an

407 argument the path to a YAML file specifying a collection of functors.

408

409 Parameters

410 ----------

411 funcs : `dict` or `list`

412 Dictionary or list of functors.

413 If a list, then it will be converted into a dictonary according to the

414 `.shortname` attribute of each functor.

415 """

416 dataset = None

417 name = "CompositeFunctor"

418

419 def __init__(self, funcs, **kwargs):

420

421 if type(funcs) is dict:

422 self.funcDict = funcs

423 else:

424 self.funcDict = {f.shortname: f for f in funcs}

425

426 self._filt = None

427

428 super().__init__(**kwargs)

429

430 @property

431 def filt(self):

432 return self._filt

433

434 @filt.setter

435 def filt(self, filt):

436 if filt is not None:

437 for _, f in self.funcDict.items():

438 f.filt = filt

439 self._filt = filt

440

441 def update(self, new):

442 """Update the functor with new functors."""

443 if isinstance(new, dict):

444 self.funcDict.update(new)

445 elif isinstance(new, CompositeFunctor):

446 self.funcDict.update(new.funcDict)

447 else:

448 raise TypeError('Can only update with dictionary or CompositeFunctor.')

449

450 # Make sure new functors have the same 'filt' set.

451 if self.filt is not None:

452 self.filt = self.filt

453

454 @property

455 def columns(self):

456 return list(set([x for y in [f.columns for f in self.funcDict.values()] for x in y]))

457

458 def multilevelColumns(self, data, **kwargs):

459 # Get the union of columns for all component functors.

460 # Note the need to have `returnTuple=True` here.

461 return list(

462 set(

463 [

464 x

465 for y in [

466 f.multilevelColumns(data, returnTuple=True, **kwargs) for f in self.funcDict.values()

467 ]

468 for x in y

469 ]

470 )

471 )

472

473 def __call__(self, data, **kwargs):

474 """Apply the functor to the data table.

475

476 Parameters

477 ----------

478 data : various

479 The data represented as `~lsst.daf.butler.DeferredDatasetHandle`,

480 `~lsst.pipe.base.InMemoryDatasetHandle`, or `~pandas.DataFrame`.

481 The table or a pointer to a table on disk from which columns can

482 be accessed.

483 """

484 if isinstance(data, pd.DataFrame):

485 _data = InMemoryDatasetHandle(data, storageClass="DataFrame")

486 elif isinstance(data, (DeferredDatasetHandle, InMemoryDatasetHandle)):

487 _data = data

488 else:

489 raise RuntimeError(f"Unexpected type provided for data. Got {get_full_type_name(data)}.")

490

491 columnIndex = self._get_columnIndex(_data)

492

493 if isinstance(columnIndex, pd.MultiIndex):

494 columns = self.multilevelColumns(_data, columnIndex=columnIndex)

495 df = _data.get(parameters={"columns": columns})

496

497 valDict = {}

498 for k, f in self.funcDict.items():

499 try:

500 subdf = f._setLevels(

501 df[f.multilevelColumns(_data, returnTuple=True, columnIndex=columnIndex)]

502 )

503 valDict[k] = f._func(subdf)

504 except Exception as e:

505 self.log.exception(

506 "Exception in %s (funcs: %s) call: %s",

507 self.name,

508 str(list(self.funcDict.keys())),

509 type(e).__name__,

510 )

511 try:

512 valDict[k] = f.fail(subdf)

513 except NameError:

514 raise e

515

516 else:

517 df = _data.get(parameters={"columns": self.columns})

518

519 valDict = {k: f._func(df) for k, f in self.funcDict.items()}

520

521 # Check that output columns are actually columns.

522 for name, colVal in valDict.items():

523 if len(colVal.shape) != 1:

524 raise RuntimeError("Transformed column '%s' is not the shape of a column. "

525 "It is shaped %s and type %s." % (name, colVal.shape, type(colVal)))

526

527 try:

528 valDf = pd.concat(valDict, axis=1)

529 except TypeError:

530 print([(k, type(v)) for k, v in valDict.items()])

531 raise

532

533 if kwargs.get('dropna', False):

534 valDf = valDf.dropna(how='any')

535

536 return valDf

537

538 @classmethod

539 def renameCol(cls, col, renameRules):

540 if renameRules is None:

541 return col

542 for old, new in renameRules:

543 if col.startswith(old):

544 col = col.replace(old, new)

545 return col

546

547 @classmethod

548 def from_file(cls, filename, **kwargs):

549 # Allow environment variables in the filename.

550 filename = os.path.expandvars(filename)

551 with open(filename) as f:

552 translationDefinition = yaml.safe_load(f)

553

554 return cls.from_yaml(translationDefinition, **kwargs)

555

556 @classmethod

557 def from_yaml(cls, translationDefinition, **kwargs):

558 funcs = {}

559 for func, val in translationDefinition['funcs'].items():

560 funcs[func] = init_fromDict(val, name=func)

561

562 if 'flag_rename_rules' in translationDefinition:

563 renameRules = translationDefinition['flag_rename_rules']

564 else:

565 renameRules = None

566

567 if 'calexpFlags' in translationDefinition:

568 for flag in translationDefinition['calexpFlags']:

569 funcs[cls.renameCol(flag, renameRules)] = Column(flag, dataset='calexp')

570

571 if 'refFlags' in translationDefinition:

572 for flag in translationDefinition['refFlags']:

573 funcs[cls.renameCol(flag, renameRules)] = Column(flag, dataset='ref')

574

575 if 'forcedFlags' in translationDefinition:

576 for flag in translationDefinition['forcedFlags']:

577 funcs[cls.renameCol(flag, renameRules)] = Column(flag, dataset='forced_src')

578

579 if 'flags' in translationDefinition:

580 for flag in translationDefinition['flags']:

581 funcs[cls.renameCol(flag, renameRules)] = Column(flag, dataset='meas')

582

583 return cls(funcs, **kwargs)

584

585

586def mag_aware_eval(df, expr, log):

587 """Evaluate an expression on a DataFrame, knowing what the 'mag' function

588 means.

589

590 Builds on `pandas.DataFrame.eval`, which parses and executes math on

591 DataFrames.

592

593 Parameters

594 ----------

595 df : ~pandas.DataFrame

596 DataFrame on which to evaluate expression.

597

598 expr : str

599 Expression.

600 """

601 try:

602 expr_new = re.sub(r'mag$(\w+)$', r'-2.5*log(\g<1>)/log(10)', expr)

603 val = df.eval(expr_new)

604 except Exception as e: # Should check what actually gets raised

605 log.error("Exception in mag_aware_eval: %s: %s", type(e).__name__, e)

606 expr_new = re.sub(r'mag$(\w+)$', r'-2.5*log(\g<1>_instFlux)/log(10)', expr)

607 val = df.eval(expr_new)

608 return val

609

610

611class CustomFunctor(Functor):

612 """Arbitrary computation on a catalog.

613

614 Column names (and thus the columns to be loaded from catalog) are found by

615 finding all words and trying to ignore all "math-y" words.

616

617 Parameters

618 ----------

619 expr : str

620 Expression to evaluate, to be parsed and executed by

621 `~lsst.pipe.tasks.functors.mag_aware_eval`.

622 """

623 _ignore_words = ('mag', 'sin', 'cos', 'exp', 'log', 'sqrt')

624

625 def __init__(self, expr, **kwargs):

626 self.expr = expr

627 super().__init__(**kwargs)

628

629 @property

630 def name(self):

631 return self.expr

632

633 @property

634 def columns(self):

635 flux_cols = re.findall(r'mag$\s*(\w+)\s*$', self.expr)

636

637 cols = [c for c in re.findall(r'[a-zA-Z_]+', self.expr) if c not in self._ignore_words]

638 not_a_col = []

639 for c in flux_cols:

640 if not re.search('_instFlux$', c):

641 cols.append(f'{c}_instFlux')

642 not_a_col.append(c)

643 else:

644 cols.append(c)

645

646 return list(set([c for c in cols if c not in not_a_col]))

647

648 def _func(self, df):

649 return mag_aware_eval(df, self.expr, self.log)

650

651

652class Column(Functor):

653 """Get column with a specified name."""

654

655 def __init__(self, col, **kwargs):

656 self.col = col

657 super().__init__(**kwargs)

658

659 @property

660 def name(self):

661 return self.col

662

663 @property

664 def columns(self):

665 return [self.col]

666

667 def _func(self, df):

668 return df[self.col]

669

670

671class Index(Functor):

672 """Return the value of the index for each object."""

673

674 columns = ['coord_ra'] # Just a dummy; something has to be here.

675 _defaultDataset = 'ref'

676 _defaultNoDup = True

677

678 def _func(self, df):

679 return pd.Series(df.index, index=df.index)

680

681

682class CoordColumn(Column):

683 """Base class for coordinate column, in degrees."""

684 _radians = True

685

686 def __init__(self, col, **kwargs):

687 super().__init__(col, **kwargs)

688

689 def _func(self, df):

690 # Must not modify original column in case that column is used by

691 # another functor.

692 output = df[self.col] * 180 / np.pi if self._radians else df[self.col]

693 return output

694

695

696class RAColumn(CoordColumn):

697 """Right Ascension, in degrees."""

698 name = 'RA'

699 _defaultNoDup = True

700

701 def __init__(self, **kwargs):

702 super().__init__('coord_ra', **kwargs)

703

704 def __call__(self, catalog, **kwargs):

705 return super().__call__(catalog, **kwargs)

706

707

708class DecColumn(CoordColumn):

709 """Declination, in degrees."""

710 name = 'Dec'

711 _defaultNoDup = True

712

713 def __init__(self, **kwargs):

714 super().__init__('coord_dec', **kwargs)

715

716 def __call__(self, catalog, **kwargs):

717 return super().__call__(catalog, **kwargs)

718

719

720class RAErrColumn(CoordColumn):

721 """Angular uncertainty in Right Ascension, in degrees.

722

723 Note that ``coord_raErr`` is the uncertainty on RA·cos(Dec),

724 i.e. an arc-length on the sky in the RA direction.

725 """

726 name = 'RAErr'

727 _defaultNoDup = True

728

729 def __init__(self, **kwargs):

730 super().__init__('coord_raErr', **kwargs)

731

732

733class DecErrColumn(CoordColumn):

734 """Uncertainty in declination, in degrees."""

735 name = 'DecErr'

736 _defaultNoDup = True

737

738 def __init__(self, **kwargs):

739 super().__init__('coord_decErr', **kwargs)

740

741

742class RADecCovColumn(Column):

743 """Tangent-plane angular covariance, in degrees^2.

744

745 As with the RA error, this is the covariance between RA·cos(Dec) and Dec,

746 Cov(xi, eta) = Cov(RA*cos(Dec), Dec).

747 """

748 _radians = True

749 name = 'RADecCov'

750 _defaultNoDup = True

751

752 def __init__(self, **kwargs):

753 super().__init__('coord_ra_dec_Cov', **kwargs)

754

755 def _func(self, df):

756 # Must not modify original column in case that column is used by

757 # another functor.

758 output = df[self.col]*(180/np.pi)**2 if self._radians else df[self.col]

759 return output

760

761

762class MultibandColumn(Column):

763 """A column with a band in a multiband table."""

764 def __init__(self, col, band_to_check, **kwargs):

765 self._band_to_check = band_to_check

766 super().__init__(col=col, **kwargs)

767

768 @property

769 def band_to_check(self):

770 return self._band_to_check

771

772

773class MultibandSinglePrecisionFloatColumn(MultibandColumn):

774 """A float32 MultibandColumn"""

775 def _func(self, df):

776 return super()._func(df).astype(np.float32)

777

778

779class SinglePrecisionFloatColumn(Column):

780 """Return a column cast to a single-precision float."""

781

782 def _func(self, df):

783 return df[self.col].astype(np.float32)

784

785

786class HtmIndex20(Functor):

787 """Compute the level 20 HtmIndex for the catalog.

788

789 Notes

790 -----

791 This functor was implemented to satisfy requirements of old APDB interface

792 which required the ``pixelId`` column in DiaObject with HTM20 index.

793 The APDB interface had migrated to not need that information, but we keep

794 this class in case it may be useful for something else.

795 """

796 name = "Htm20"

797 htmLevel = 20

798 _radians = True

799

800 def __init__(self, ra, dec, **kwargs):

801 self.pixelator = sphgeom.HtmPixelization(self.htmLevel)

802 self.ra = ra

803 self.dec = dec

804 self._columns = [self.ra, self.dec]

805 super().__init__(**kwargs)

806

807 def _func(self, df):

808

809 def computePixel(row):

810 if self._radians:

811 sphPoint = geom.SpherePoint(row[self.ra],

812 row[self.dec],

813 geom.radians)

814 else:

815 sphPoint = geom.SpherePoint(row[self.ra],

816 row[self.dec],

817 geom.degrees)

818 return self.pixelator.index(sphPoint.getVector())

819

820 return df.apply(computePixel, axis=1, result_type='reduce').astype('int64')

821

822

823def fluxName(col):

824 """Append _instFlux to the column name if it doesn't have it already."""

825 if not col.endswith('_instFlux'):

826 col += '_instFlux'

827 return col

828

829

830def fluxErrName(col):

831 """Append _instFluxErr to the column name if it doesn't have it already."""

832 if not col.endswith('_instFluxErr'):

833 col += '_instFluxErr'

834 return col

835

836

837class Mag(Functor):

838 """Compute calibrated magnitude.

839

840 Returns the flux at mag=0.

841 The default ``fluxMag0`` is 63095734448.0194, which is default for HSC.

842 TO DO: This default should be made configurable in DM-21955.

843

844 This calculation hides warnings about invalid values and dividing by zero.

845

846 As with all functors, a ``dataset`` and ``filt`` kwarg should be provided

847 upon initialization.

848 Unlike the default `Functor`, however, the default dataset for a `Mag` is

849 ``'meas'``, rather than ``'ref'``.

850

851 Parameters

852 ----------

853 col : `str`

854 Name of flux column from which to compute magnitude.

855 Can be parseable by the `~lsst.pipe.tasks.functors.fluxName` function;

856 that is, you can pass ``'modelfit_CModel'`` instead of

857 ``'modelfit_CModel_instFlux'``, and it will understand.

858 """

859 _defaultDataset = 'meas'

860

861 def __init__(self, col, **kwargs):

862 self.col = fluxName(col)

863 # TO DO: DM-21955 Replace hard coded photometic calibration values.

864 self.fluxMag0 = 63095734448.0194

865

866 super().__init__(**kwargs)

867

868 @property

869 def columns(self):

870 return [self.col]

871

872 def _func(self, df):

873 with warnings.catch_warnings():

874 warnings.filterwarnings('ignore', r'invalid value encountered')

875 warnings.filterwarnings('ignore', r'divide by zero')

876 return -2.5*np.log10(df[self.col] / self.fluxMag0)

877

878 @property

879 def name(self):

880 return f'mag_{self.col}'

881

882

883class MagErr(Mag):

884 """Compute calibrated magnitude uncertainty.

885

886 Parameters

887 ----------

888 col : `str`

889 Name of the flux column.

890 """

891

892 def __init__(self, *args, **kwargs):

893 super().__init__(*args, **kwargs)

894 # TO DO: DM-21955 Replace hard coded photometic calibration values.

895 self.fluxMag0Err = 0.

896

897 @property

898 def columns(self):

899 return [self.col, self.col + 'Err']

900

901 def _func(self, df):

902 with warnings.catch_warnings():

903 warnings.filterwarnings('ignore', r'invalid value encountered')

904 warnings.filterwarnings('ignore', r'divide by zero')

905 fluxCol, fluxErrCol = self.columns

906 x = df[fluxErrCol] / df[fluxCol]

907 y = self.fluxMag0Err / self.fluxMag0

908 magErr = (2.5 / np.log(10.)) * np.sqrt(x*x + y*y)

909 return magErr

910

911 @property

912 def name(self):

913 return super().name + '_err'

914

915

916class MagDiff(Functor):

917 """Functor to calculate magnitude difference."""

918 _defaultDataset = 'meas'

919

920 def __init__(self, col1, col2, **kwargs):

921 self.col1 = fluxName(col1)

922 self.col2 = fluxName(col2)

923 super().__init__(**kwargs)

924

925 @property

926 def columns(self):

927 return [self.col1, self.col2]

928

929 def _func(self, df):

930 with warnings.catch_warnings():

931 warnings.filterwarnings('ignore', r'invalid value encountered')

932 warnings.filterwarnings('ignore', r'divide by zero')

933 return -2.5*np.log10(df[self.col1]/df[self.col2])

934

935 @property

936 def name(self):

937 return f'(mag_{self.col1} - mag_{self.col2})'

938

939 @property

940 def shortname(self):

941 return f'magDiff_{self.col1}_{self.col2}'

942

943

944class Color(Functor):

945 """Compute the color between two filters.

946

947 Computes color by initializing two different `Mag` functors based on the

948 ``col`` and filters provided, and then returning the difference.

949

950 This is enabled by the `_func` method expecting a DataFrame with a

951 multilevel column index, with both ``'band'`` and ``'column'``, instead of

952 just ``'column'``, which is the `Functor` default.

953 This is controlled by the `_dfLevels` attribute.

954

955 Also of note, the default dataset for `Color` is ``forced_src'``, whereas

956 for `Mag` it is ``'meas'``.

957

958 Parameters

959 ----------

960 col : str

961 Name of the flux column from which to compute; same as would be passed

962 to `~lsst.pipe.tasks.functors.Mag`.

963

964 filt2, filt1 : str

965 Filters from which to compute magnitude difference.

966 Color computed is ``Mag(filt2) - Mag(filt1)``.

967 """

968 _defaultDataset = 'forced_src'

969 _dfLevels = ('band', 'column')

970 _defaultNoDup = True

971

972 def __init__(self, col, filt2, filt1, **kwargs):

973 self.col = fluxName(col)

974 if filt2 == filt1:

975 raise RuntimeError("Cannot compute Color for %s: %s - %s " % (col, filt2, filt1))

976 self.filt2 = filt2

977 self.filt1 = filt1

978

979 self.mag2 = Mag(col, filt=filt2, **kwargs)

980 self.mag1 = Mag(col, filt=filt1, **kwargs)

981

982 super().__init__(**kwargs)

983

984 @property

985 def filt(self):

986 return None

987

988 @filt.setter

989 def filt(self, filt):

990 pass

991

992 def _func(self, df):

993 mag2 = self.mag2._func(df[self.filt2])

994 mag1 = self.mag1._func(df[self.filt1])

995 return mag2 - mag1

996

997 @property

998 def columns(self):

999 return [self.mag1.col, self.mag2.col]

1000

1001 def multilevelColumns(self, parq, **kwargs):

1002 return [(self.dataset, self.filt1, self.col), (self.dataset, self.filt2, self.col)]

1003

1004 @property

1005 def name(self):

1006 return f'{self.filt2} - {self.filt1} ({self.col})'

1007

1008 @property

1009 def shortname(self):

1010 return f"{self.col}_{self.filt2.replace('-', '')}m{self.filt1.replace('-', '')}"

1011

1012

1013class DeconvolvedMoments(Functor):

1014 """This functor subtracts the trace of the PSF second moments from the

1015 trace of the second moments of the source.

1016

1017 If the HsmShapeAlgorithm measurement is valid, then these will be used for

1018 the sources.

1019 Otherwise, the SdssShapeAlgorithm measurements will be used.

1020 """

1021 name = 'Deconvolved Moments'

1022 shortname = 'deconvolvedMoments'

1023 _columns = ("ext_shapeHSM_HsmSourceMoments_xx",

1024 "ext_shapeHSM_HsmSourceMoments_yy",

1025 "base_SdssShape_xx", "base_SdssShape_yy",

1026 "ext_shapeHSM_HsmPsfMoments_xx",

1027 "ext_shapeHSM_HsmPsfMoments_yy")

1028

1029 def _func(self, df):

1030 """Calculate deconvolved moments."""

1031 if "ext_shapeHSM_HsmSourceMoments_xx" in df.columns: # _xx added by tdm

1032 hsm = df["ext_shapeHSM_HsmSourceMoments_xx"] + df["ext_shapeHSM_HsmSourceMoments_yy"]

1033 else:

1034 hsm = np.ones(len(df))*np.nan

1035 sdss = df["base_SdssShape_xx"] + df["base_SdssShape_yy"]

1036 if "ext_shapeHSM_HsmPsfMoments_xx" in df.columns:

1037 psf = df["ext_shapeHSM_HsmPsfMoments_xx"] + df["ext_shapeHSM_HsmPsfMoments_yy"]

1038 else:

1039 # LSST does not have shape.sdss.psf.

1040 # We could instead add base_PsfShape to the catalog using

1041 # exposure.getPsf().computeShape(s.getCentroid()).getIxx().

1042 raise RuntimeError('No psf shape parameter found in catalog')

1043

1044 return hsm.where(np.isfinite(hsm), sdss) - psf

1045

1046

1047class SdssTraceSize(Functor):

1048 """Functor to calculate the SDSS trace radius size for sources.

1049

1050 The SDSS trace radius size is a measure of size equal to the square root of

1051 half of the trace of the second moments tensor measured with the

1052 SdssShapeAlgorithm plugin.

1053 This has units of pixels.

1054 """

1055 name = "SDSS Trace Size"

1056 shortname = 'sdssTrace'

1057 _columns = ("base_SdssShape_xx", "base_SdssShape_yy")

1058

1059 def _func(self, df):

1060 srcSize = np.sqrt(0.5*(df["base_SdssShape_xx"] + df["base_SdssShape_yy"]))

1061 return srcSize

1062

1063

1064class PsfSdssTraceSizeDiff(Functor):

1065 """Functor to calculate the SDSS trace radius size difference (%) between

1066 the object and the PSF model.

1067

1068 See Also

1069 --------

1070 SdssTraceSize

1071 """

1072 name = "PSF - SDSS Trace Size"

1073 shortname = 'psf_sdssTrace'

1074 _columns = ("base_SdssShape_xx", "base_SdssShape_yy",

1075 "base_SdssShape_psf_xx", "base_SdssShape_psf_yy")

1076

1077 def _func(self, df):

1078 srcSize = np.sqrt(0.5*(df["base_SdssShape_xx"] + df["base_SdssShape_yy"]))

1079 psfSize = np.sqrt(0.5*(df["base_SdssShape_psf_xx"] + df["base_SdssShape_psf_yy"]))

1080 sizeDiff = 100*(srcSize - psfSize)/(0.5*(srcSize + psfSize))

1081 return sizeDiff

1082

1083

1084class HsmTraceSize(Functor):

1085 """Functor to calculate the HSM trace radius size for sources.

1086

1087 The HSM trace radius size is a measure of size equal to the square root of

1088 half of the trace of the second moments tensor measured with the

1089 HsmShapeAlgorithm plugin.

1090 This has units of pixels.

1091 """

1092 name = 'HSM Trace Size'

1093 shortname = 'hsmTrace'

1094 _columns = ("ext_shapeHSM_HsmSourceMoments_xx",

1095 "ext_shapeHSM_HsmSourceMoments_yy")

1096

1097 def _func(self, df):

1098 srcSize = np.sqrt(0.5*(df["ext_shapeHSM_HsmSourceMoments_xx"]

1099 + df["ext_shapeHSM_HsmSourceMoments_yy"]))

1100 return srcSize

1101

1102

1103class PsfHsmTraceSizeDiff(Functor):

1104 """Functor to calculate the HSM trace radius size difference (%) between

1105 the object and the PSF model.

1106

1107 See Also

1108 --------

1109 HsmTraceSize

1110 """

1111 name = 'PSF - HSM Trace Size'

1112 shortname = 'psf_HsmTrace'

1113 _columns = ("ext_shapeHSM_HsmSourceMoments_xx",

1114 "ext_shapeHSM_HsmSourceMoments_yy",

1115 "ext_shapeHSM_HsmPsfMoments_xx",

1116 "ext_shapeHSM_HsmPsfMoments_yy")

1117

1118 def _func(self, df):

1119 srcSize = np.sqrt(0.5*(df["ext_shapeHSM_HsmSourceMoments_xx"]

1120 + df["ext_shapeHSM_HsmSourceMoments_yy"]))

1121 psfSize = np.sqrt(0.5*(df["ext_shapeHSM_HsmPsfMoments_xx"]

1122 + df["ext_shapeHSM_HsmPsfMoments_yy"]))

1123 sizeDiff = 100*(srcSize - psfSize)/(0.5*(srcSize + psfSize))

1124 return sizeDiff

1125

1126

1127class HsmFwhm(Functor):

1128 """Functor to calculate the PSF FWHM with second moments measured from the

1129 HsmShapeAlgorithm plugin.

1130

1131 This is in units of arcseconds, and assumes the hsc_rings_v1 skymap pixel

1132 scale of 0.168 arcseconds/pixel.

1133

1134 Notes

1135 -----

1136 This conversion assumes the PSF is Gaussian, which is not always the case.

1137 """

1138 name = 'HSM Psf FWHM'

1139 _columns = ('ext_shapeHSM_HsmPsfMoments_xx', 'ext_shapeHSM_HsmPsfMoments_yy')

1140 # TODO: DM-21403 pixel scale should be computed from the CD matrix or transform matrix

1141 pixelScale = 0.168

1142 SIGMA2FWHM = 2*np.sqrt(2*np.log(2))

1143

1144 def _func(self, df):

1145 return (self.pixelScale*self.SIGMA2FWHM*np.sqrt(

1146 0.5*(df['ext_shapeHSM_HsmPsfMoments_xx']

1147 + df['ext_shapeHSM_HsmPsfMoments_yy']))).astype(np.float32)

1148

1149

1150class E1(Functor):

1151 r"""Calculate :math:`e_1` ellipticity component for sources, defined as:

1152

1153 .. math::

1154 e_1 &= (I_{xx}-I_{yy})/(I_{xx}+I_{yy})

1155

1156 See Also

1157 --------

1158 E2

1159 """

1160 name = "Distortion Ellipticity (e1)"

1161 shortname = "Distortion"

1162

1163 def __init__(self, colXX, colXY, colYY, **kwargs):

1164 self.colXX = colXX

1165 self.colXY = colXY

1166 self.colYY = colYY

1167 self._columns = [self.colXX, self.colXY, self.colYY]

1168 super().__init__(**kwargs)

1169

1170 @property

1171 def columns(self):

1172 return [self.colXX, self.colXY, self.colYY]

1173

1174 def _func(self, df):

1175 return ((df[self.colXX] - df[self.colYY]) / (

1176 df[self.colXX] + df[self.colYY])).astype(np.float32)

1177

1178

1179class E2(Functor):

1180 r"""Calculate :math:`e_2` ellipticity component for sources, defined as:

1181

1182 .. math::

1183 e_2 &= 2I_{xy}/(I_{xx}+I_{yy})

1184

1185 See Also

1186 --------

1187 E1

1188 """

1189 name = "Ellipticity e2"

1190

1191 def __init__(self, colXX, colXY, colYY, **kwargs):

1192 self.colXX = colXX

1193 self.colXY = colXY

1194 self.colYY = colYY

1195 super().__init__(**kwargs)

1196

1197 @property

1198 def columns(self):

1199 return [self.colXX, self.colXY, self.colYY]

1200

1201 def _func(self, df):

1202 return (2*df[self.colXY] / (df[self.colXX] + df[self.colYY])).astype(np.float32)

1203

1204

1205class RadiusFromQuadrupole(Functor):

1206 """Calculate the radius from the quadrupole moments.

1207

1208 This returns the fourth root of the determinant of the second moments

1209 tensor, which has units of pixels.

1210

1211 See Also

1212 --------

1213 SdssTraceSize

1214 HsmTraceSize

1215 """

1216

1217 def __init__(self, colXX, colXY, colYY, **kwargs):

1218 self.colXX = colXX

1219 self.colXY = colXY

1220 self.colYY = colYY

1221 super().__init__(**kwargs)

1222

1223 @property

1224 def columns(self):

1225 return [self.colXX, self.colXY, self.colYY]

1226

1227 def _func(self, df):

1228 return ((df[self.colXX]*df[self.colYY] - df[self.colXY]**2)**0.25).astype(np.float32)

1229

1230

1231class LocalWcs(Functor):

1232 """Computations using the stored localWcs."""

1233 name = "LocalWcsOperations"

1234

1235 def __init__(self,

1236 colCD_1_1,

1237 colCD_1_2,

1238 colCD_2_1,

1239 colCD_2_2,

1240 **kwargs):

1241 self.colCD_1_1 = colCD_1_1

1242 self.colCD_1_2 = colCD_1_2

1243 self.colCD_2_1 = colCD_2_1

1244 self.colCD_2_2 = colCD_2_2

1245 super().__init__(**kwargs)

1246

1247 def computeDeltaRaDec(self, x, y, cd11, cd12, cd21, cd22):

1248 """Compute the dRA, dDec from dx, dy.

1249

1250 Parameters

1251 ----------

1252 x : `~pandas.Series`

1253 X pixel coordinate.

1254 y : `~pandas.Series`

1255 Y pixel coordinate.

1256 cd11 : `~pandas.Series`

1257 [1, 1] element of the local Wcs affine transform.

1258 cd12 : `~pandas.Series`

1259 [1, 2] element of the local Wcs affine transform.

1260 cd21 : `~pandas.Series`

1261 [2, 1] element of the local Wcs affine transform.

1262 cd22 : `~pandas.Series`

1263 [2, 2] element of the local Wcs affine transform.

1264

1265 Returns

1266 -------

1267 raDecTuple : tuple

1268 RA and Dec conversion of x and y given the local Wcs.

1269 Returned units are in radians.

1270

1271 Notes

1272 -----

1273 If x and y are with respect to the CRVAL1, CRVAL2

1274 then this will return the RA, Dec for that WCS.

1275 """

1276 return (x * cd11 + y * cd12, x * cd21 + y * cd22)

1277

1278 def computeSkySeparation(self, ra1, dec1, ra2, dec2):

1279 """Compute the local pixel scale conversion.

1280

1281 Parameters

1282 ----------

1283 ra1 : `~pandas.Series`

1284 Ra of the first coordinate in radians.

1285 dec1 : `~pandas.Series`

1286 Dec of the first coordinate in radians.

1287 ra2 : `~pandas.Series`

1288 Ra of the second coordinate in radians.

1289 dec2 : `~pandas.Series`

1290 Dec of the second coordinate in radians.

1291

1292 Returns

1293 -------

1294 dist : `~pandas.Series`

1295 Distance on the sphere in radians.

1296 """

1297 deltaDec = dec2 - dec1

1298 deltaRa = ra2 - ra1

1299 return 2 * np.arcsin(

1300 np.sqrt(

1301 np.sin(deltaDec / 2) ** 2

1302 + np.cos(dec2) * np.cos(dec1) * np.sin(deltaRa / 2) ** 2))

1303

1304 def getSkySeparationFromPixel(self, x1, y1, x2, y2, cd11, cd12, cd21, cd22):

1305 """Compute the distance on the sphere from x2, y1 to x1, y1.

1306

1307 Parameters

1308 ----------

1309 x1 : `~pandas.Series`

1310 X pixel coordinate.

1311 y1 : `~pandas.Series`

1312 Y pixel coordinate.

1313 x2 : `~pandas.Series`

1314 X pixel coordinate.

1315 y2 : `~pandas.Series`

1316 Y pixel coordinate.

1317 cd11 : `~pandas.Series`

1318 [1, 1] element of the local Wcs affine transform.

1319 cd12 : `~pandas.Series`

1320 [1, 2] element of the local Wcs affine transform.

1321 cd21 : `~pandas.Series`

1322 [2, 1] element of the local Wcs affine transform.

1323 cd22 : `~pandas.Series`

1324 [2, 2] element of the local Wcs affine transform.

1325

1326 Returns

1327 -------

1328 Distance : `~pandas.Series`

1329 Arcseconds per pixel at the location of the local WC.

1330 """

1331 ra1, dec1 = self.computeDeltaRaDec(x1, y1, cd11, cd12, cd21, cd22)

1332 ra2, dec2 = self.computeDeltaRaDec(x2, y2, cd11, cd12, cd21, cd22)

1333 # Great circle distance for small separations.

1334 return self.computeSkySeparation(ra1, dec1, ra2, dec2)

1335

1336 def computePositionAngle(self, ra1, dec1, ra2, dec2):

1337 """Compute position angle (E of N) from (ra1, dec1) to (ra2, dec2).

1338

1339 Parameters

1340 ----------

1341 ra1 : iterable [`float`]

1342 RA of the first coordinate [radian].

1343 dec1 : iterable [`float`]

1344 Dec of the first coordinate [radian].

1345 ra2 : iterable [`float`]

1346 RA of the second coordinate [radian].

1347 dec2 : iterable [`float`]

1348 Dec of the second coordinate [radian].

1349

1350 Returns

1351 -------

1352 Position Angle: `~pandas.Series`

1353 radians E of N

1354

1355 Notes

1356 -----

1357 (ra1, dec1) -> (ra2, dec2) is interpreted as the shorter way around the sphere

1358

1359 For a separation of 0.0001 rad, the position angle is good to 0.0009 rad

1360 all over the sphere.

1361 """

1362 # lsst.geom.SpherePoint has "bearingTo", which returns angle N of E

1363 # We instead want the astronomy convention of "Position Angle", which is angle E of N

1364 position_angle = np.zeros(len(ra1))

1365 for i, (r1, d1, r2, d2) in enumerate(zip(ra1, dec1, ra2, dec2)):

1366 point1 = geom.SpherePoint(r1, d1, geom.radians)

1367 point2 = geom.SpherePoint(r2, d2, geom.radians)

1368 bearing = point1.bearingTo(point2)

1369 pa_ref_angle = geom.Angle(np.pi/2, geom.radians) # in bearing system

1370 pa = pa_ref_angle - bearing

1371 # Wrap around to get Delta_RA from -pi to +pi

1372 pa = pa.wrapCtr()

1373 position_angle[i] = pa.asRadians()

1374

1375 return pd.Series(position_angle)

1376

1377 def getPositionAngleFromDetectorAngle(self, theta, cd11, cd12, cd21, cd22):

1378 """Compute position angle (E of N) from detector angle (+y of +x).

1379

1380 Parameters

1381 ----------

1382 theta : `float`

1383 detector angle [radian]

1384 cd11 : `float`

1385 [1, 1] element of the local Wcs affine transform.

1386 cd12 : `float`

1387 [1, 2] element of the local Wcs affine transform.

1388 cd21 : `float`

1389 [2, 1] element of the local Wcs affine transform.

1390 cd22 : `float`

1391 [2, 2] element of the local Wcs affine transform.

1392

1393 Returns

1394 -------

1395 Position Angle: `~pandas.Series`

1396 Degrees E of N.

1397 """

1398 # Create a unit vector in (x, y) along da

1399 dx = np.cos(theta)

1400 dy = np.sin(theta)

1401 ra1, dec1 = self.computeDeltaRaDec(0, 0, cd11, cd12, cd21, cd22)

1402 ra2, dec2 = self.computeDeltaRaDec(dx, dy, cd11, cd12, cd21, cd22)

1403 # Position angle of vector from (RA1, Dec1) to (RA2, Dec2)

1404 return np.rad2deg(self.computePositionAngle(ra1, dec1, ra2, dec2))

1405

1406

1407class ComputePixelScale(LocalWcs):

1408 """Compute the local pixel scale from the stored CDMatrix.

1409 """

1410 name = "PixelScale"

1411

1412 @property

1413 def columns(self):

1414 return [self.colCD_1_1,

1415 self.colCD_1_2,

1416 self.colCD_2_1,

1417 self.colCD_2_2]

1418

1419 def pixelScaleArcseconds(self, cd11, cd12, cd21, cd22):

1420 """Compute the local pixel to scale conversion in arcseconds.

1421

1422 Parameters

1423 ----------

1424 cd11 : `~pandas.Series`

1425 [1, 1] element of the local Wcs affine transform in radians.

1426 cd11 : `~pandas.Series`

1427 [1, 1] element of the local Wcs affine transform in radians.

1428 cd12 : `~pandas.Series`

1429 [1, 2] element of the local Wcs affine transform in radians.

1430 cd21 : `~pandas.Series`

1431 [2, 1] element of the local Wcs affine transform in radians.

1432 cd22 : `~pandas.Series`

1433 [2, 2] element of the local Wcs affine transform in radians.

1434

1435 Returns

1436 -------

1437 pixScale : `~pandas.Series`

1438 Arcseconds per pixel at the location of the local WC.

1439 """

1440 return 3600 * np.degrees(np.sqrt(np.fabs(cd11 * cd22 - cd12 * cd21)))

1441

1442 def _func(self, df):

1443 return self.pixelScaleArcseconds(df[self.colCD_1_1],

1444 df[self.colCD_1_2],

1445 df[self.colCD_2_1],

1446 df[self.colCD_2_2])

1447

1448

1449class ConvertPixelToArcseconds(ComputePixelScale):

1450 """Convert a value in units of pixels to units of arcseconds."""

1451

1452 def __init__(self,

1453 col,

1454 colCD_1_1,

1455 colCD_1_2,

1456 colCD_2_1,

1457 colCD_2_2,

1458 **kwargs):

1459 self.col = col

1460 super().__init__(colCD_1_1,

1461 colCD_1_2,

1462 colCD_2_1,

1463 colCD_2_2,

1464 **kwargs)

1465

1466 @property

1467 def name(self):

1468 return f"{self.col}_asArcseconds"

1469

1470 @property

1471 def columns(self):

1472 return [self.col,

1473 self.colCD_1_1,

1474 self.colCD_1_2,

1475 self.colCD_2_1,

1476 self.colCD_2_2]

1477

1478 def _func(self, df):

1479 return df[self.col] * self.pixelScaleArcseconds(df[self.colCD_1_1],

1480 df[self.colCD_1_2],

1481 df[self.colCD_2_1],

1482 df[self.colCD_2_2])

1483

1484

1485class ConvertPixelSqToArcsecondsSq(ComputePixelScale):

1486 """Convert a value in units of pixels squared to units of arcseconds

1487 squared.

1488 """

1489

1490 def __init__(self,

1491 col,

1492 colCD_1_1,

1493 colCD_1_2,

1494 colCD_2_1,

1495 colCD_2_2,

1496 **kwargs):

1497 self.col = col

1498 super().__init__(colCD_1_1,

1499 colCD_1_2,

1500 colCD_2_1,

1501 colCD_2_2,

1502 **kwargs)

1503

1504 @property

1505 def name(self):

1506 return f"{self.col}_asArcsecondsSq"

1507

1508 @property

1509 def columns(self):

1510 return [self.col,

1511 self.colCD_1_1,

1512 self.colCD_1_2,

1513 self.colCD_2_1,

1514 self.colCD_2_2]

1515

1516 def _func(self, df):

1517 pixScale = self.pixelScaleArcseconds(df[self.colCD_1_1],

1518 df[self.colCD_1_2],

1519 df[self.colCD_2_1],

1520 df[self.colCD_2_2])

1521 return df[self.col] * pixScale * pixScale

1522

1523

1524class ConvertDetectorAngleToPositionAngle(LocalWcs):

1525 """Compute a position angle from a detector angle and the stored CDMatrix.

1526

1527 Returns

1528 -------

1529 position angle : degrees

1530 """

1531

1532 name = "PositionAngle"

1533

1534 def __init__(

1535 self,

1536 theta_col,

1537 colCD_1_1,

1538 colCD_1_2,

1539 colCD_2_1,

1540 colCD_2_2,

1541 **kwargs

1542 ):

1543 self.theta_col = theta_col

1544 super().__init__(colCD_1_1, colCD_1_2, colCD_2_1, colCD_2_2, **kwargs)

1545

1546 @property

1547 def columns(self):

1548 return [

1549 self.theta_col,

1550 self.colCD_1_1,

1551 self.colCD_1_2,

1552 self.colCD_2_1,

1553 self.colCD_2_2

1554 ]

1555

1556 def _func(self, df):

1557 return self.getPositionAngleFromDetectorAngle(

1558 df[self.theta_col],

1559 df[self.colCD_1_1],

1560 df[self.colCD_1_2],

1561 df[self.colCD_2_1],

1562 df[self.colCD_2_2]

1563 )

1564

1565

1566class ConvertDetectorAngleErrToPositionAngleErr(Functor):

1567 """Convert a detector angle error to a position angle error.

1568

1569 Returns

1570 -------

1571 position angle error : degrees

1572 """

1573

1574 name = "PositionAngleErr"

1575

1576 def __init__(self, theta_err_col, **kwargs):

1577 self.theta_err_col = theta_err_col

1578 super().__init__(**kwargs)

1579

1580 @property

1581 def columns(self):

1582 return [self.theta_err_col]

1583

1584 def _func(self, df):

1585 return np.rad2deg(df[self.theta_err_col])

1586

1587

1588class ReferenceBand(Functor):

1589 """Return the band used to seed multiband forced photometry.

1590

1591 This functor is to be used on Object tables.

1592 It converts the boolean merge_measurements_{band} columns into a single

1593 string representing the first band for which merge_measurements_{band}

1594 is True.

1595

1596 Assumes the default priority order of i, r, z, y, g, u.

1597 """

1598 name = 'Reference Band'

1599 shortname = 'refBand'

1600

1601 band_order = ("i", "r", "z", "y", "g", "u")

1602

1603 @property

1604 def columns(self):

1605 # Build the actual input column list, not hardcoded ugrizy

1606 bands = [band for band in self.band_order if band in self.bands]

1607 # In the unlikely scenario that users attempt to add non-ugrizy bands

1608 bands += [band for band in self.bands if band not in self.band_order]

1609 return [f"merge_measurement_{band}" for band in bands]

1610

1611 def _func(self, df: pd.DataFrame) -> pd.Series:

1612 def getFilterAliasName(row):

1613 # Get column name with the max value (True > False).

1614 colName = row.idxmax()

1615 return colName.replace('merge_measurement_', '')

1616

1617 # Skip columns that are unavailable, because this functor requests the

1618 # superset of bands that could be included in the object table.

1619 columns = [col for col in self.columns if col in df.columns]

1620 # Makes a Series of dtype object if df is empty.

1621 return df[columns].apply(getFilterAliasName, axis=1,

1622 result_type='reduce').astype('object')

1623

1624 def __init__(self, bands: tuple[str] | list[str] | None = None, **kwargs):

1625 super().__init__(**kwargs)

1626 self.bands = self.band_order if bands is None else tuple(bands)

1627

1628

1629class Photometry(Functor):

1630 """Base class for Object table calibrated fluxes and magnitudes."""

1631 # AB to NanoJansky (3631 Jansky).

1632 AB_FLUX_SCALE = (0 * u.ABmag).to_value(u.nJy)

1633 LOG_AB_FLUX_SCALE = 12.56

1634 FIVE_OVER_2LOG10 = 1.085736204758129569

1635 # TO DO: DM-21955 Replace hard coded photometic calibration values.

1636 COADD_ZP = 27

1637

1638 def __init__(self, colFlux, colFluxErr=None, **kwargs):

1639 self.vhypot = np.vectorize(self.hypot)

1640 self.col = colFlux

1641 self.colFluxErr = colFluxErr

1642

1643 self.fluxMag0 = 1./np.power(10, -0.4*self.COADD_ZP)

1644 self.fluxMag0Err = 0.

1645

1646 super().__init__(**kwargs)

1647

1648 @property

1649 def columns(self):

1650 return [self.col]

1651

1652 @property

1653 def name(self):

1654 return f'mag_{self.col}'

1655

1656 @classmethod

1657 def hypot(cls, a, b):

1658 """Compute sqrt(a^2 + b^2) without under/overflow."""

1659 if np.abs(a) < np.abs(b):

1660 a, b = b, a

1661 if a == 0.:

1662 return 0.

1663 q = b/a

1664 return np.abs(a) * np.sqrt(1. + q*q)

1665

1666 def dn2flux(self, dn, fluxMag0):

1667 """Convert instrumental flux to nanojanskys."""

1668 return (self.AB_FLUX_SCALE * dn / fluxMag0).astype(np.float32)

1669

1670 def dn2mag(self, dn, fluxMag0):

1671 """Convert instrumental flux to AB magnitude."""

1672 with warnings.catch_warnings():

1673 warnings.filterwarnings('ignore', r'invalid value encountered')

1674 warnings.filterwarnings('ignore', r'divide by zero')

1675 return (-2.5 * np.log10(dn/fluxMag0)).astype(np.float32)

1676

1677 def dn2fluxErr(self, dn, dnErr, fluxMag0, fluxMag0Err):

1678 """Convert instrumental flux error to nanojanskys."""

1679 retVal = self.vhypot(dn * fluxMag0Err, dnErr * fluxMag0)

1680 retVal *= self.AB_FLUX_SCALE / fluxMag0 / fluxMag0

1681 return retVal.astype(np.float32)

1682

1683 def dn2MagErr(self, dn, dnErr, fluxMag0, fluxMag0Err):

1684 """Convert instrumental flux error to AB magnitude error."""

1685 retVal = self.dn2fluxErr(dn, dnErr, fluxMag0, fluxMag0Err) / self.dn2flux(dn, fluxMag0)

1686 return (self.FIVE_OVER_2LOG10 * retVal).astype(np.float32)

1687

1688

1689class NanoJansky(Photometry):

1690 """Convert instrumental flux to nanojanskys."""

1691 def _func(self, df):

1692 return self.dn2flux(df[self.col], self.fluxMag0)

1693

1694

1695class NanoJanskyErr(Photometry):

1696 """Convert instrumental flux error to nanojanskys."""

1697 @property

1698 def columns(self):

1699 return [self.col, self.colFluxErr]

1700

1701 def _func(self, df):

1702 retArr = self.dn2fluxErr(df[self.col], df[self.colFluxErr], self.fluxMag0, self.fluxMag0Err)

1703 return pd.Series(retArr, index=df.index)

1704

1705

1706class LocalPhotometry(Functor):

1707 """Base class for calibrating the specified instrument flux column using

1708 the local photometric calibration.

1709

1710 Parameters

1711 ----------

1712 instFluxCol : `str`

1713 Name of the instrument flux column.

1714 instFluxErrCol : `str`

1715 Name of the assocated error columns for ``instFluxCol``.

1716 photoCalibCol : `str`

1717 Name of local calibration column.

1718 photoCalibErrCol : `str`, optional

1719 Error associated with ``photoCalibCol``. Ignored and deprecated; will

1720 be removed after v29.

1721

1722 See Also

1723 --------

1724 LocalNanojansky

1725 LocalNanojanskyErr

1726 """

1727 logNJanskyToAB = (1 * u.nJy).to_value(u.ABmag)

1728

1729 def __init__(self,

1730 instFluxCol,

1731 instFluxErrCol,

1732 photoCalibCol,

1733 photoCalibErrCol=None,

1734 **kwargs):

1735 self.instFluxCol = instFluxCol

1736 self.instFluxErrCol = instFluxErrCol

1737 self.photoCalibCol = photoCalibCol

1738 # TODO[DM-49400]: remove this check and the argument it corresponds to.

1739 if photoCalibErrCol is not None:

1740 warnings.warn("The photoCalibErrCol argument is deprecated and will be removed after v29.",

1741 category=FutureWarning)

1742 super().__init__(**kwargs)

1743

1744 def instFluxToNanojansky(self, instFlux, localCalib):

1745 """Convert instrument flux to nanojanskys.

1746

1747 Parameters

1748 ----------

1749 instFlux : `~numpy.ndarray` or `~pandas.Series`

1750 Array of instrument flux measurements.

1751 localCalib : `~numpy.ndarray` or `~pandas.Series`

1752 Array of local photometric calibration estimates.

1753

1754 Returns

1755 -------

1756 calibFlux : `~numpy.ndarray` or `~pandas.Series`

1757 Array of calibrated flux measurements.

1758 """

1759 return instFlux * localCalib

1760

1761 def instFluxErrToNanojanskyErr(self, instFlux, instFluxErr, localCalib, localCalibErr=None):

1762 """Convert instrument flux to nanojanskys.

1763

1764 Parameters

1765 ----------

1766 instFlux : `~numpy.ndarray` or `~pandas.Series`

1767 Array of instrument flux measurements. Ignored (accepted for

1768 backwards compatibility and consistency with magnitude-error

1769 calculation methods).

1770 instFluxErr : `~numpy.ndarray` or `~pandas.Series`

1771 Errors on associated ``instFlux`` values.

1772 localCalib : `~numpy.ndarray` or `~pandas.Series`

1773 Array of local photometric calibration estimates.

1774 localCalibErr : `~numpy.ndarray` or `~pandas.Series`, optional

1775 Errors on associated ``localCalib`` values. Ignored and deprecated;

1776 will be removed after v29.

1777

1778 Returns

1779 -------

1780 calibFluxErr : `~numpy.ndarray` or `~pandas.Series`

1781 Errors on calibrated flux measurements.

1782 """

1783 # TODO[DM-49400]: remove this check and the argument it corresponds to.

1784 if localCalibErr is not None:

1785 warnings.warn("The localCalibErr argument is deprecated and will be removed after v29.",

1786 category=FutureWarning)

1787 return instFluxErr * localCalib

1788

1789 def instFluxToMagnitude(self, instFlux, localCalib):

1790 """Convert instrument flux to nanojanskys.

1791

1792 Parameters

1793 ----------

1794 instFlux : `~numpy.ndarray` or `~pandas.Series`

1795 Array of instrument flux measurements.

1796 localCalib : `~numpy.ndarray` or `~pandas.Series`

1797 Array of local photometric calibration estimates.

1798

1799 Returns

1800 -------

1801 calibMag : `~numpy.ndarray` or `~pandas.Series`

1802 Array of calibrated AB magnitudes.

1803 """

1804 return -2.5 * np.log10(self.instFluxToNanojansky(instFlux, localCalib)) + self.logNJanskyToAB

1805

1806 def instFluxErrToMagnitudeErr(self, instFlux, instFluxErr, localCalib, localCalibErr=None):

1807 """Convert instrument flux err to nanojanskys.

1808

1809 Parameters

1810 ----------

1811 instFlux : `~numpy.ndarray` or `~pandas.Series`

1812 Array of instrument flux measurements.

1813 instFluxErr : `~numpy.ndarray` or `~pandas.Series`

1814 Errors on associated ``instFlux`` values.

1815 localCalib : `~numpy.ndarray` or `~pandas.Series`

1816 Array of local photometric calibration estimates.

1817 localCalibErr : `~numpy.ndarray` or `~pandas.Series`, optional

1818 Errors on associated ``localCalib`` values. Ignored and deprecated;

1819 will be removed after v29.

1820

1821 Returns

1822 -------

1823 calibMagErr: `~numpy.ndarray` or `~pandas.Series`

1824 Error on calibrated AB magnitudes.

1825 """

1826 # TODO[DM-49400]: remove this check and the argument it corresponds to.

1827 if localCalibErr is not None:

1828 warnings.warn("The localCalibErr argument is deprecated and will be removed after v29.",

1829 category=FutureWarning)

1830 err = self.instFluxErrToNanojanskyErr(instFlux, instFluxErr, localCalib)

1831 return 2.5 / np.log(10) * err / self.instFluxToNanojansky(instFlux, instFluxErr)

1832

1833

1834class LocalNanojansky(LocalPhotometry):

1835 """Compute calibrated fluxes using the local calibration value.

1836

1837 This returns units of nanojanskys.

1838 """

1839

1840 @property

1841 def columns(self):

1842 return [self.instFluxCol, self.photoCalibCol]

1843

1844 @property

1845 def name(self):

1846 return f'flux_{self.instFluxCol}'

1847

1848 def _func(self, df):

1849 return self.instFluxToNanojansky(df[self.instFluxCol],

1850 df[self.photoCalibCol]).astype(np.float32)

1851

1852

1853class LocalNanojanskyErr(LocalPhotometry):

1854 """Compute calibrated flux errors using the local calibration value.

1855

1856 This returns units of nanojanskys.

1857 """

1858

1859 @property

1860 def columns(self):

1861 return [self.instFluxCol, self.instFluxErrCol, self.photoCalibCol]

1862

1863 @property

1864 def name(self):

1865 return f'fluxErr_{self.instFluxCol}'

1866

1867 def _func(self, df):

1868 return self.instFluxErrToNanojanskyErr(df[self.instFluxCol], df[self.instFluxErrCol],

1869 df[self.photoCalibCol]).astype(np.float32)

1870

1871

1872class LocalDipoleMeanFlux(LocalPhotometry):

1873 """Compute absolute mean of dipole fluxes.

1874

1875 See Also

1876 --------

1877 LocalNanojansky

1878 LocalNanojanskyErr

1879 LocalDipoleMeanFluxErr

1880 LocalDipoleDiffFlux

1881 LocalDipoleDiffFluxErr

1882 """

1883 def __init__(self,

1884 instFluxPosCol,

1885 instFluxNegCol,

1886 instFluxPosErrCol,

1887 instFluxNegErrCol,

1888 photoCalibCol,

1889 # TODO[DM-49400]: remove this option; it's already deprecated (in super).

1890 photoCalibErrCol=None,

1891 **kwargs):

1892 self.instFluxNegCol = instFluxNegCol

1893 self.instFluxPosCol = instFluxPosCol

1894 self.instFluxNegErrCol = instFluxNegErrCol

1895 self.instFluxPosErrCol = instFluxPosErrCol

1896 self.photoCalibCol = photoCalibCol

1897 super().__init__(instFluxNegCol,

1898 instFluxNegErrCol,

1899 photoCalibCol,

1900 photoCalibErrCol,

1901 **kwargs)

1902

1903 @property

1904 def columns(self):

1905 return [self.instFluxPosCol,

1906 self.instFluxNegCol,

1907 self.photoCalibCol]

1908

1909 @property

1910 def name(self):

1911 return f'dipMeanFlux_{self.instFluxPosCol}_{self.instFluxNegCol}'

1912

1913 def _func(self, df):

1914 return 0.5*(np.fabs(self.instFluxToNanojansky(df[self.instFluxNegCol], df[self.photoCalibCol]))

1915 + np.fabs(self.instFluxToNanojansky(df[self.instFluxPosCol], df[self.photoCalibCol])))

1916

1917

1918class LocalDipoleMeanFluxErr(LocalDipoleMeanFlux):

1919 """Compute the error on the absolute mean of dipole fluxes.

1920

1921 See Also

1922 --------

1923 LocalNanojansky

1924 LocalNanojanskyErr

1925 LocalDipoleMeanFlux

1926 LocalDipoleDiffFlux

1927 LocalDipoleDiffFluxErr

1928 """

1929

1930 @property

1931 def columns(self):

1932 return [self.instFluxPosCol,

1933 self.instFluxNegCol,

1934 self.instFluxPosErrCol,

1935 self.instFluxNegErrCol,

1936 self.photoCalibCol]

1937

1938 @property

1939 def name(self):

1940 return f'dipMeanFluxErr_{self.instFluxPosCol}_{self.instFluxNegCol}'

1941

1942 def _func(self, df):

1943 return 0.5*np.hypot(df[self.instFluxNegErrCol], df[self.instFluxPosErrCol]) * df[self.photoCalibCol]

1944

1945

1946class LocalDipoleDiffFlux(LocalDipoleMeanFlux):

1947 """Compute the absolute difference of dipole fluxes.

1948

1949 Calculated value is (abs(pos) - abs(neg)).

1950

1951 See Also

1952 --------

1953 LocalNanojansky

1954 LocalNanojanskyErr

1955 LocalDipoleMeanFlux

1956 LocalDipoleMeanFluxErr

1957 LocalDipoleDiffFluxErr

1958 """

1959

1960 @property

1961 def columns(self):

1962 return [self.instFluxPosCol,

1963 self.instFluxNegCol,

1964 self.photoCalibCol]

1965

1966 @property

1967 def name(self):

1968 return f'dipDiffFlux_{self.instFluxPosCol}_{self.instFluxNegCol}'

1969

1970 def _func(self, df):

1971 return (np.fabs(self.instFluxToNanojansky(df[self.instFluxPosCol], df[self.photoCalibCol]))

1972 - np.fabs(self.instFluxToNanojansky(df[self.instFluxNegCol], df[self.photoCalibCol])))

1973

1974

1975class LocalDipoleDiffFluxErr(LocalDipoleMeanFlux):

1976 """Compute the error on the absolute difference of dipole fluxes.

1977

1978 See Also

1979 --------

1980 LocalNanojansky

1981 LocalNanojanskyErr

1982 LocalDipoleMeanFlux

1983 LocalDipoleMeanFluxErr

1984 LocalDipoleDiffFlux

1985 """

1986

1987 @property

1988 def columns(self):

1989 return [self.instFluxPosCol,

1990 self.instFluxNegCol,

1991 self.instFluxPosErrCol,

1992 self.instFluxNegErrCol,

1993 self.photoCalibCol]

1994

1995 @property

1996 def name(self):

1997 return f'dipDiffFluxErr_{self.instFluxPosCol}_{self.instFluxNegCol}'

1998

1999 def _func(self, df):

2000 return np.hypot(df[self.instFluxPosErrCol], df[self.instFluxNegErrCol]) * df[self.photoCalibCol]

2001

2002

2003class Ebv(Functor):

2004 """Compute E(B-V) from dustmaps.sfd."""

2005 _defaultDataset = 'ref'

2006 name = "E(B-V)"

2007 shortname = "ebv"

2008

2009 def __init__(self, **kwargs):

2010 # Import is only needed for Ebv.

2011 # Suppress unnecessary .dustmapsrc log message on import.

2012 with open(os.devnull, "w") as devnull:

2013 with redirect_stdout(devnull):

2014 from dustmaps.sfd import SFDQuery

2015 self._columns = ['coord_ra', 'coord_dec']

2016 self.sfd = SFDQuery()

2017 super().__init__(**kwargs)

2018

2019 def _func(self, df):

2020 coords = SkyCoord(df['coord_ra'].values * u.rad, df['coord_dec'].values * u.rad)

2021 ebv = self.sfd(coords)

2022 return pd.Series(ebv, index=df.index).astype('float32')

2023

2024

2025class MomentsBase(Functor):

2026 """Base class for functors that use shape moments and localWCS

2027

2028 Attributes

2029 ----------

2030 is_covariance : bool

2031 Whether the shape columns are terms of a covariance matrix. If False,

2032 they will be assumed to be terms of a correlation matrix instead.

2033 """

2034

2035 is_covariance: bool = True

2036

2037 def __init__(self,

2038 shape_1_1,

2039 shape_2_2,

2040 shape_1_2,

2041 colCD_1_1,

2042 colCD_1_2,

2043 colCD_2_1,

2044 colCD_2_2,

2045 **kwargs):

2046 self.shape_1_1 = shape_1_1

2047 self.shape_2_2 = shape_2_2

2048 self.shape_1_2 = shape_1_2

2049 self.colCD_1_1 = colCD_1_1

2050 self.colCD_1_2 = colCD_1_2

2051 self.colCD_2_1 = colCD_2_1

2052 self.colCD_2_2 = colCD_2_2

2053 super().__init__(**kwargs)

2054

2055 @property

2056 def columns(self):

2057 return [

2058 self.shape_1_1,

2059 self.shape_2_2,

2060 self.shape_1_2,

2061 ] + self.columns_ref

2062

2063 @property

2064 def columns_ref(self):

2065 """Return columns that are needed from the ref table."""

2066 return [

2067 self.colCD_1_1,

2068 self.colCD_1_2,

2069 self.colCD_2_1,

2070 self.colCD_2_2]

2071

2072 def compute_ellipse_terms(self, df, sky: bool = True):

2073 r"""Return terms commonly used for ellipse parameterization conversions.

2074

2075 Parameters

2076 ----------

2077 df

2078 The data frame.

2079 sky

2080 Whether to compute the terms in sky coordinates.

2081 If False, XX, YY and XY moments are used instead of

2082 UU, VV and UV.

2083

2084 Returns

2085 -------

2086 xx_p_yy

2087 The sum of the diagonal terms of the covariance.

2088 xx_m_yy

2089 The difference of the diagonal terms of the covariance.

2090 t2

2091 A term similar to the discriminant of the quadratic formula.

2092 """

2093 xx = self.sky_uu(df) if sky else self.get_xx(df)

2094 yy = self.sky_vv(df) if sky else self.get_yy(df)

2095 xx_m_yy = xx - yy

2096 t2 = xx_m_yy**2 + 4.0*(self.sky_uv(df) if sky else self.get_xy(df))**2

2097 # TODO: Check alternative form that may be more stable for computing

2098 # the minor axis size (see gauss2d/src/ellipse.cc)

2099 # t2 = xx**2 + yy**2 - 2*(xx*yy - 2*xy**2)

2100 return xx + yy, xx_m_yy, t2

2101

2102 def get_xx(self, df):

2103 xx = df[self.shape_1_1]

2104 return xx if self.is_covariance else xx**2

2105

2106 def get_yy(self, df):

2107 yy = df[self.shape_2_2]

2108 return yy if self.is_covariance else yy**2

2109

2110 def get_xy(self, df):

2111 xy = df[self.shape_1_2]

2112 return xy if self.is_covariance else xy*df[self.shape_1_1]*df[self.shape_2_2]

2113

2114 # Each of sky_uu, sky_vv, sky_uv evaluates one element of

2115 # CD_matrix * moments_matrix * CD_matrix.T

2116 def sky_uu(self, df):

2117 """Return the component of the moments tensor aligned with the RA axis, in radians."""

2118 i_xx = self.get_xx(df)

2119 i_yy = self.get_yy(df)

2120 i_xy = self.get_xy(df)

2121 CD_1_1 = df[self.colCD_1_1]

2122 CD_1_2 = df[self.colCD_1_2]

2123 return (CD_1_1*(i_xx*CD_1_1 + i_xy*CD_1_2)

2124 + CD_1_2*(i_xy*CD_1_1 + i_yy*CD_1_2))

2125

2126 def sky_vv(self, df):

2127 """Return the component of the moments tensor aligned with the dec axis, in radians."""

2128 i_xx = self.get_xx(df)

2129 i_yy = self.get_yy(df)

2130 i_xy = self.get_xy(df)

2131 CD_2_1 = df[self.colCD_2_1]

2132 CD_2_2 = df[self.colCD_2_2]

2133 return (CD_2_1*(i_xx*CD_2_1 + i_xy*CD_2_2)

2134 + CD_2_2*(i_xy*CD_2_1 + i_yy*CD_2_2))

2135

2136 def sky_uv(self, df):

2137 """Return the covariance of the moments tensor in ra, dec coordinates, in radians."""

2138 i_xx = self.get_xx(df)

2139 i_yy = self.get_yy(df)

2140 i_xy = self.get_xy(df)

2141 CD_1_1 = df[self.colCD_1_1]

2142 CD_1_2 = df[self.colCD_1_2]

2143 CD_2_1 = df[self.colCD_2_1]

2144 CD_2_2 = df[self.colCD_2_2]

2145 return ((CD_1_1 * i_xx + CD_1_2 * i_xy) * CD_2_1

2146 + (CD_1_1 * i_xy + CD_1_2 * i_yy) * CD_2_2)

2147

2148 def get_g1(self, df):

2149 """

2150 Calculate shear-type ellipticity parameter G1.

2151 """

2152 # TODO: Replace this with functionality from afwGeom, DM-54015

2153 sky_uu = self.sky_uu(df)

2154 sky_vv = self.sky_vv(df)

2155 sky_uv = self.sky_uv(df)

2156 denom = sky_uu + sky_vv + 2 * np.sqrt(sky_uu*sky_vv - sky_uv**2)

2157 return ((sky_uu - sky_vv) / denom).astype(np.float32)

2158

2159 def get_g2(self, df):

2160 """

2161 Calculate shear-type ellipticity parameter G2.

2162

2163 This has the opposite sign as sky_uv in order to maintain consistency with the HSM moments

2164 sign convention.

2165 """

2166 # TODO: Replace this with functionality from afwGeom, DM-54015

2167 sky_uu = self.sky_uu(df)

2168 sky_vv = self.sky_vv(df)

2169 sky_uv = self.sky_uv(df)

2170 denom = sky_uu + sky_vv + 2 * np.sqrt(sky_uu*sky_vv - sky_uv**2)

2171 return (-2*sky_uv / denom).astype(np.float32)

2172

2173 def get_trace(self, df):

2174 sky_uu = self.sky_uu(df)

2175 sky_vv = self.sky_vv(df)

2176 return np.sqrt(0.5*(sky_uu + sky_vv)).astype(np.float32)

2177

2178

2179class MomentsG1Sky(MomentsBase):

2180 """Rotate pixel moments Ixx,Iyy,Ixy into RA/dec frame and G1/G2 reduced

2181 shear parameterization"""

2182 _defaultDataset = 'meas'

2183 name = "moments_g1"

2184 shortname = "moments_g1"

2185

2186 def _func(self, df):

2187 sky_g1 = self.get_g1(df)

2188

2189 return pd.Series(sky_g1.astype(np.float32), index=df.index)

2190

2191

2192class MomentsG2Sky(MomentsBase):

2193 """Rotate pixel moments Ixx,Iyy,Ixy into RA/dec frame and G1/G2 reduced

2194 shear parameterization"""

2195 _defaultDataset = 'meas'

2196 name = "moments_g2"

2197 shortname = "moments_g2"

2198

2199 def _func(self, df):

2200 sky_g2 = self.get_g2(df)

2201

2202 return pd.Series(sky_g2.astype(np.float32), index=df.index)

2203

2204

2205class MomentsTraceSky(MomentsBase):

2206 """Trace radius size in arcseconds from pixel moments Ixx,Iyy,Ixy

2207

2208 The trace radius size is a measure of size equal to the square root of

2209 half of the trace of the second moments tensor.

2210 """

2211 _defaultDataset = 'meas'

2212 name = "moments_trace"

2213 shortname = "moments_trace"

2214

2215 def _func(self, df):

2216 sky_trace_radians = self.get_trace(df)

2217

2218 return pd.Series((sky_trace_radians*(180/np.pi)*3600).astype(np.float32), index=df.index)

2219

2220

2221class MomentsIuuSky(MomentsBase):

2222 """Rotate pixel moments Ixx,Iyy,Ixy into ra,dec frame and arcseconds"""

2223 _defaultDataset = 'meas'

2224 name = "moments_uu"

2225 shortname = "moments_uu"

2226

2227 def _func(self, df):

2228 sky_uu_radians = self.sky_uu(df)

2229

2230 return pd.Series((sky_uu_radians*((180/np.pi)*3600)**2).astype(np.float32), index=df.index)

2231

2232

2233class CorrelationIuuSky(MomentsIuuSky):

2234 """MomentsIuuSky but from sigma_x, sigma_y, rho correlation terms."""

2235 is_covariance = False

2236

2237

2238class MomentsIvvSky(MomentsBase):

2239 """Rotate pixel moments Ixx,Iyy,Ixy into ra,dec frame and arcseconds"""

2240 _defaultDataset = 'meas'

2241 name = "moments_vv"

2242 shortname = "moments_vv"

2243

2244 def _func(self, df):

2245 sky_vv_radians = self.sky_vv(df)

2246

2247 return pd.Series((sky_vv_radians*((180/np.pi)*3600)**2).astype(np.float32), index=df.index)

2248

2249

2250class CorrelationIvvSky(MomentsIvvSky):

2251 """MomentsIvvSky but from sigma_x, sigma_y, rho correlation terms."""

2252 is_covariance = False

2253

2254

2255class MomentsIuvSky(MomentsBase):

2256 """Rotate pixel moments Ixx,Iyy,Ixy into ra,dec frame and arcseconds"""

2257 _defaultDataset = 'meas'

2258 name = "moments_uv"

2259 shortname = "moments_uv"

2260

2261 def _func(self, df):

2262 sky_uv_radians = self.sky_uv(df)

2263

2264 return pd.Series((sky_uv_radians*((180/np.pi)*3600)**2).astype(np.float32), index=df.index)

2265

2266

2267class CorrelationIuvSky(MomentsIuvSky):

2268 """MomentsIuvSky but from sigma_x, sigma_y, rho correlation terms."""

2269 is_covariance = False

2270

2271

2272class PositionAngleFromMoments(MomentsBase):

2273 """Compute position angle relative to ra,dec frame, in degrees, from Ixx,Iyy,Ixy pixel moments."""

2274 _defaultDataset = 'meas'

2275 name = "moments_theta"

2276 shortname = "moments_theta"

2277

2278 def _func(self, df):

2279 sky_uu = self.sky_uu(df)

2280 sky_vv = self.sky_vv(df)

2281 sky_uv = self.sky_uv(df)

2282 theta = 0.5*np.arctan2(2*sky_uv, sky_uu - sky_vv)

2283

2284 return pd.Series((np.degrees(np.array(theta))).astype(np.float32), index=df.index)

2285

2286

2287class PositionAngleFromCorrelation(PositionAngleFromMoments):

2288 """PositionAngleFromMoments but from sigma_x, sigma_y, rho correlation terms."""

2289 is_covariance = False

2290

2291

2292class SemimajorAxisFromMoments(MomentsBase):

2293 """Compute the semimajor axis length in arcseconds, from Ixx,Iyy,Ixy pixel moments."""

2294 _defaultDataset = 'meas'

2295 name = "moments_a"

2296 shortname = "moments_a"

2297

2298 def _func(self, df):

2299 xx_p_yy, _, t2 = self.compute_ellipse_terms(df)

2300 # This copies what is done (unvectorized) in afw.geom.ellipse

2301 a_radians = np.sqrt(0.5 * (xx_p_yy + np.sqrt(t2)))

2302

2303 return pd.Series((np.degrees(a_radians)*3600).astype(np.float32), index=df.index)

2304

2305

2306class SemimajorAxisFromCorrelation(SemimajorAxisFromMoments):

2307 """SemimajorAxisFromMoments but from sigma_x, sigma_y, rho correlation terms."""

2308 is_covariance = False

2309

2310

2311class SemiminorAxisFromMoments(MomentsBase):

2312 """Compute the semiminor axis length in arcseconds, from Ixx,Iyy,Ixy pixel moments."""

2313 _defaultDataset = 'meas'

2314 name = "moments_b"

2315 shortname = "moments_b"

2316

2317 def _func(self, df):

2318 xx_p_yy, _, t2 = self.compute_ellipse_terms(df)

2319 # This copies what is done (unvectorized) in afw.geom.ellipse

2320 b_radians = np.sqrt(0.5 * (xx_p_yy - np.sqrt(t2)))

2321

2322 return pd.Series((np.degrees(b_radians)*3600).astype(np.float32), index=df.index)

2323

2324

2325class SemiminorAxisFromCorrelation(SemiminorAxisFromMoments):

2326 """SemiminorAxisFromMoments but from sigma_x, sigma_y, rho correlation terms."""

2327 is_covariance = False