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 """Uncertainty in Right Ascension, in degrees."""

722 name = 'RAErr'

723 _defaultNoDup = True

724

725 def __init__(self, **kwargs):

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

727

728

729class DecErrColumn(CoordColumn):

730 """Uncertainty in declination, in degrees."""

731 name = 'DecErr'

732 _defaultNoDup = True

733

734 def __init__(self, **kwargs):

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

736

737

738class RADecCovColumn(Column):

739 """Coordinate covariance column, in degrees."""

740 _radians = True

741 name = 'RADecCov'

742 _defaultNoDup = True

743

744 def __init__(self, **kwargs):

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

746

747 def _func(self, df):

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

749 # another functor.

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

751 return output

752

753

754class MultibandColumn(Column):

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

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

757 self._band_to_check = band_to_check

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

759

760 @property

761 def band_to_check(self):

762 return self._band_to_check

763

764

765class MultibandSinglePrecisionFloatColumn(MultibandColumn):

766 """A float32 MultibandColumn"""

767 def _func(self, df):

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

769

770

771class SinglePrecisionFloatColumn(Column):

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

773

774 def _func(self, df):

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

776

777

778class HtmIndex20(Functor):

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

780

781 Notes

782 -----

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

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

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

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

787 """

788 name = "Htm20"

789 htmLevel = 20

790 _radians = True

791

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

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

794 self.ra = ra

795 self.dec = dec

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

797 super().__init__(**kwargs)

798

799 def _func(self, df):

800

801 def computePixel(row):

802 if self._radians:

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

804 row[self.dec],

805 geom.radians)

806 else:

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

808 row[self.dec],

809 geom.degrees)

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

811

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

813

814

815def fluxName(col):

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

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

818 col += '_instFlux'

819 return col

820

821

822def fluxErrName(col):

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

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

825 col += '_instFluxErr'

826 return col

827

828

829class Mag(Functor):

830 """Compute calibrated magnitude.

831

832 Returns the flux at mag=0.

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

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

835

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

837

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

839 upon initialization.

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

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

842

843 Parameters

844 ----------

845 col : `str`

846 Name of flux column from which to compute magnitude.

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

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

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

850 """

851 _defaultDataset = 'meas'

852

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

854 self.col = fluxName(col)

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

856 self.fluxMag0 = 63095734448.0194

857

858 super().__init__(**kwargs)

859

860 @property

861 def columns(self):

862 return [self.col]

863

864 def _func(self, df):

865 with warnings.catch_warnings():

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

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

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

869

870 @property

871 def name(self):

872 return f'mag_{self.col}'

873

874

875class MagErr(Mag):

876 """Compute calibrated magnitude uncertainty.

877

878 Parameters

879 ----------

880 col : `str`

881 Name of the flux column.

882 """

883

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

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

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

887 self.fluxMag0Err = 0.

888

889 @property

890 def columns(self):

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

892

893 def _func(self, df):

894 with warnings.catch_warnings():

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

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

897 fluxCol, fluxErrCol = self.columns

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

899 y = self.fluxMag0Err / self.fluxMag0

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

901 return magErr

902

903 @property

904 def name(self):

905 return super().name + '_err'

906

907

908class MagDiff(Functor):

909 """Functor to calculate magnitude difference."""

910 _defaultDataset = 'meas'

911

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

913 self.col1 = fluxName(col1)

914 self.col2 = fluxName(col2)

915 super().__init__(**kwargs)

916

917 @property

918 def columns(self):

919 return [self.col1, self.col2]

920

921 def _func(self, df):

922 with warnings.catch_warnings():

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

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

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

926

927 @property

928 def name(self):

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

930

931 @property

932 def shortname(self):

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

934

935

936class Color(Functor):

937 """Compute the color between two filters.

938

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

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

941

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

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

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

945 This is controlled by the `_dfLevels` attribute.

946

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

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

949

950 Parameters

951 ----------

952 col : str

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

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

955

956 filt2, filt1 : str

957 Filters from which to compute magnitude difference.

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

959 """

960 _defaultDataset = 'forced_src'

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

962 _defaultNoDup = True

963

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

965 self.col = fluxName(col)

966 if filt2 == filt1:

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

968 self.filt2 = filt2

969 self.filt1 = filt1

970

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

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

973

974 super().__init__(**kwargs)

975

976 @property

977 def filt(self):

978 return None

979

980 @filt.setter

981 def filt(self, filt):

982 pass

983

984 def _func(self, df):

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

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

987 return mag2 - mag1

988

989 @property

990 def columns(self):

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

992

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

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

995

996 @property

997 def name(self):

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

999

1000 @property

1001 def shortname(self):

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

1003

1004

1005class DeconvolvedMoments(Functor):

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

1007 trace of the second moments of the source.

1008

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

1010 the sources.

1011 Otherwise, the SdssShapeAlgorithm measurements will be used.

1012 """

1013 name = 'Deconvolved Moments'

1014 shortname = 'deconvolvedMoments'

1015 _columns = ("ext_shapeHSM_HsmSourceMoments_xx",

1016 "ext_shapeHSM_HsmSourceMoments_yy",

1017 "base_SdssShape_xx", "base_SdssShape_yy",

1018 "ext_shapeHSM_HsmPsfMoments_xx",

1019 "ext_shapeHSM_HsmPsfMoments_yy")

1020

1021 def _func(self, df):

1022 """Calculate deconvolved moments."""

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

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

1025 else:

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

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

1028 if "ext_shapeHSM_HsmPsfMoments_xx" in df.columns:

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

1030 else:

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

1032 # We could instead add base_PsfShape to the catalog using

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

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

1035

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

1037

1038

1039class SdssTraceSize(Functor):

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

1041

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

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

1044 SdssShapeAlgorithm plugin.

1045 This has units of pixels.

1046 """

1047 name = "SDSS Trace Size"

1048 shortname = 'sdssTrace'

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

1050

1051 def _func(self, df):

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

1053 return srcSize

1054

1055

1056class PsfSdssTraceSizeDiff(Functor):

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

1058 the object and the PSF model.

1059

1060 See Also

1061 --------

1062 SdssTraceSize

1063 """

1064 name = "PSF - SDSS Trace Size"

1065 shortname = 'psf_sdssTrace'

1066 _columns = ("base_SdssShape_xx", "base_SdssShape_yy",

1067 "base_SdssShape_psf_xx", "base_SdssShape_psf_yy")

1068

1069 def _func(self, df):

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

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

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

1073 return sizeDiff

1074

1075

1076class HsmTraceSize(Functor):

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

1078

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

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

1081 HsmShapeAlgorithm plugin.

1082 This has units of pixels.

1083 """

1084 name = 'HSM Trace Size'

1085 shortname = 'hsmTrace'

1086 _columns = ("ext_shapeHSM_HsmSourceMoments_xx",

1087 "ext_shapeHSM_HsmSourceMoments_yy")

1088

1089 def _func(self, df):

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

1091 + df["ext_shapeHSM_HsmSourceMoments_yy"]))

1092 return srcSize

1093

1094

1095class PsfHsmTraceSizeDiff(Functor):

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

1097 the object and the PSF model.

1098

1099 See Also

1100 --------

1101 HsmTraceSize

1102 """

1103 name = 'PSF - HSM Trace Size'

1104 shortname = 'psf_HsmTrace'

1105 _columns = ("ext_shapeHSM_HsmSourceMoments_xx",

1106 "ext_shapeHSM_HsmSourceMoments_yy",

1107 "ext_shapeHSM_HsmPsfMoments_xx",

1108 "ext_shapeHSM_HsmPsfMoments_yy")

1109

1110 def _func(self, df):

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

1112 + df["ext_shapeHSM_HsmSourceMoments_yy"]))

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

1114 + df["ext_shapeHSM_HsmPsfMoments_yy"]))

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

1116 return sizeDiff

1117

1118

1119class HsmFwhm(Functor):

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

1121 HsmShapeAlgorithm plugin.

1122

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

1124 scale of 0.168 arcseconds/pixel.

1125

1126 Notes

1127 -----

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

1129 """

1130 name = 'HSM Psf FWHM'

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

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

1133 pixelScale = 0.168

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

1135

1136 def _func(self, df):

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

1138 0.5*(df['ext_shapeHSM_HsmPsfMoments_xx']

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

1140

1141

1142class E1(Functor):

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

1144

1145 .. math::

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

1147

1148 See Also

1149 --------

1150 E2

1151 """

1152 name = "Distortion Ellipticity (e1)"

1153 shortname = "Distortion"

1154

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

1156 self.colXX = colXX

1157 self.colXY = colXY

1158 self.colYY = colYY

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

1160 super().__init__(**kwargs)

1161

1162 @property

1163 def columns(self):

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

1165

1166 def _func(self, df):

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

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

1169

1170

1171class E2(Functor):

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

1173

1174 .. math::

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

1176

1177 See Also

1178 --------

1179 E1

1180 """

1181 name = "Ellipticity e2"

1182

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

1184 self.colXX = colXX

1185 self.colXY = colXY

1186 self.colYY = colYY

1187 super().__init__(**kwargs)

1188

1189 @property

1190 def columns(self):

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

1192

1193 def _func(self, df):

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

1195

1196

1197class RadiusFromQuadrupole(Functor):

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

1199

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

1201 tensor, which has units of pixels.

1202

1203 See Also

1204 --------

1205 SdssTraceSize

1206 HsmTraceSize

1207 """

1208

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

1210 self.colXX = colXX

1211 self.colXY = colXY

1212 self.colYY = colYY

1213 super().__init__(**kwargs)

1214

1215 @property

1216 def columns(self):

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

1218

1219 def _func(self, df):

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

1221

1222

1223class LocalWcs(Functor):

1224 """Computations using the stored localWcs."""

1225 name = "LocalWcsOperations"

1226

1227 def __init__(self,

1228 colCD_1_1,

1229 colCD_1_2,

1230 colCD_2_1,

1231 colCD_2_2,

1232 **kwargs):

1233 self.colCD_1_1 = colCD_1_1

1234 self.colCD_1_2 = colCD_1_2

1235 self.colCD_2_1 = colCD_2_1

1236 self.colCD_2_2 = colCD_2_2

1237 super().__init__(**kwargs)

1238

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

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

1241

1242 Parameters

1243 ----------

1244 x : `~pandas.Series`

1245 X pixel coordinate.

1246 y : `~pandas.Series`

1247 Y pixel coordinate.

1248 cd11 : `~pandas.Series`

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

1250 cd12 : `~pandas.Series`

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

1252 cd21 : `~pandas.Series`

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

1254 cd22 : `~pandas.Series`

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

1256

1257 Returns

1258 -------

1259 raDecTuple : tuple

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

1261 Returned units are in radians.

1262

1263 Notes

1264 -----

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

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

1267 """

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

1269

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

1271 """Compute the local pixel scale conversion.

1272

1273 Parameters

1274 ----------

1275 ra1 : `~pandas.Series`

1276 Ra of the first coordinate in radians.

1277 dec1 : `~pandas.Series`

1278 Dec of the first coordinate in radians.

1279 ra2 : `~pandas.Series`

1280 Ra of the second coordinate in radians.

1281 dec2 : `~pandas.Series`

1282 Dec of the second coordinate in radians.

1283

1284 Returns

1285 -------

1286 dist : `~pandas.Series`

1287 Distance on the sphere in radians.

1288 """

1289 deltaDec = dec2 - dec1

1290 deltaRa = ra2 - ra1

1291 return 2 * np.arcsin(

1292 np.sqrt(

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

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

1295

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

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

1298

1299 Parameters

1300 ----------

1301 x1 : `~pandas.Series`

1302 X pixel coordinate.

1303 y1 : `~pandas.Series`

1304 Y pixel coordinate.

1305 x2 : `~pandas.Series`

1306 X pixel coordinate.

1307 y2 : `~pandas.Series`

1308 Y pixel coordinate.

1309 cd11 : `~pandas.Series`

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

1311 cd12 : `~pandas.Series`

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

1313 cd21 : `~pandas.Series`

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

1315 cd22 : `~pandas.Series`

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

1317

1318 Returns

1319 -------

1320 Distance : `~pandas.Series`

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

1322 """

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

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

1325 # Great circle distance for small separations.

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

1327

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

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

1330

1331 Parameters

1332 ----------

1333 ra1 : iterable [`float`]

1334 RA of the first coordinate [radian].

1335 dec1 : iterable [`float`]

1336 Dec of the first coordinate [radian].

1337 ra2 : iterable [`float`]

1338 RA of the second coordinate [radian].

1339 dec2 : iterable [`float`]

1340 Dec of the second coordinate [radian].

1341

1342 Returns

1343 -------

1344 Position Angle: `~pandas.Series`

1345 radians E of N

1346

1347 Notes

1348 -----

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

1350

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

1352 all over the sphere.

1353 """

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

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

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

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

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

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

1360 bearing = point1.bearingTo(point2)

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

1362 pa = pa_ref_angle - bearing

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

1364 pa = pa.wrapCtr()

1365 position_angle[i] = pa.asRadians()

1366

1367 return pd.Series(position_angle)

1368

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

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

1371

1372 Parameters

1373 ----------

1374 theta : `float`

1375 detector angle [radian]

1376 cd11 : `float`

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

1378 cd12 : `float`

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

1380 cd21 : `float`

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

1382 cd22 : `float`

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

1384

1385 Returns

1386 -------

1387 Position Angle: `~pandas.Series`

1388 Degrees E of N.

1389 """

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

1391 dx = np.cos(theta)

1392 dy = np.sin(theta)

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

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

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

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

1397

1398

1399class ComputePixelScale(LocalWcs):

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

1401 """

1402 name = "PixelScale"

1403

1404 @property

1405 def columns(self):

1406 return [self.colCD_1_1,

1407 self.colCD_1_2,

1408 self.colCD_2_1,

1409 self.colCD_2_2]

1410

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

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

1413

1414 Parameters

1415 ----------

1416 cd11 : `~pandas.Series`

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

1418 cd11 : `~pandas.Series`

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

1420 cd12 : `~pandas.Series`

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

1422 cd21 : `~pandas.Series`

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

1424 cd22 : `~pandas.Series`

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

1426

1427 Returns

1428 -------

1429 pixScale : `~pandas.Series`

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

1431 """

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

1433

1434 def _func(self, df):

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

1436 df[self.colCD_1_2],

1437 df[self.colCD_2_1],

1438 df[self.colCD_2_2])

1439

1440

1441class ConvertPixelToArcseconds(ComputePixelScale):

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

1443

1444 def __init__(self,

1445 col,

1446 colCD_1_1,

1447 colCD_1_2,

1448 colCD_2_1,

1449 colCD_2_2,

1450 **kwargs):

1451 self.col = col

1452 super().__init__(colCD_1_1,

1453 colCD_1_2,

1454 colCD_2_1,

1455 colCD_2_2,

1456 **kwargs)

1457

1458 @property

1459 def name(self):

1460 return f"{self.col}_asArcseconds"

1461

1462 @property

1463 def columns(self):

1464 return [self.col,

1465 self.colCD_1_1,

1466 self.colCD_1_2,

1467 self.colCD_2_1,

1468 self.colCD_2_2]

1469

1470 def _func(self, df):

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

1472 df[self.colCD_1_2],

1473 df[self.colCD_2_1],

1474 df[self.colCD_2_2])

1475

1476

1477class ConvertPixelSqToArcsecondsSq(ComputePixelScale):

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

1479 squared.

1480 """

1481

1482 def __init__(self,

1483 col,

1484 colCD_1_1,

1485 colCD_1_2,

1486 colCD_2_1,

1487 colCD_2_2,

1488 **kwargs):

1489 self.col = col

1490 super().__init__(colCD_1_1,

1491 colCD_1_2,

1492 colCD_2_1,

1493 colCD_2_2,

1494 **kwargs)

1495

1496 @property

1497 def name(self):

1498 return f"{self.col}_asArcsecondsSq"

1499

1500 @property

1501 def columns(self):

1502 return [self.col,

1503 self.colCD_1_1,

1504 self.colCD_1_2,

1505 self.colCD_2_1,

1506 self.colCD_2_2]

1507

1508 def _func(self, df):

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

1510 df[self.colCD_1_2],

1511 df[self.colCD_2_1],

1512 df[self.colCD_2_2])

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

1514

1515

1516class ConvertDetectorAngleToPositionAngle(LocalWcs):

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

1518

1519 Returns

1520 -------

1521 position angle : degrees

1522 """

1523

1524 name = "PositionAngle"

1525

1526 def __init__(

1527 self,

1528 theta_col,

1529 colCD_1_1,

1530 colCD_1_2,

1531 colCD_2_1,

1532 colCD_2_2,

1533 **kwargs

1534 ):

1535 self.theta_col = theta_col

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

1537

1538 @property

1539 def columns(self):

1540 return [

1541 self.theta_col,

1542 self.colCD_1_1,

1543 self.colCD_1_2,

1544 self.colCD_2_1,

1545 self.colCD_2_2

1546 ]

1547

1548 def _func(self, df):

1549 return self.getPositionAngleFromDetectorAngle(

1550 df[self.theta_col],

1551 df[self.colCD_1_1],

1552 df[self.colCD_1_2],

1553 df[self.colCD_2_1],

1554 df[self.colCD_2_2]

1555 )

1556

1557

1558class ConvertDetectorAngleErrToPositionAngleErr(Functor):

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

1560

1561 Returns

1562 -------

1563 position angle error : degrees

1564 """

1565

1566 name = "PositionAngleErr"

1567

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

1569 self.theta_err_col = theta_err_col

1570 super().__init__(**kwargs)

1571

1572 @property

1573 def columns(self):

1574 return [self.theta_err_col]

1575

1576 def _func(self, df):

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

1578

1579

1580class ReferenceBand(Functor):

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

1582

1583 This functor is to be used on Object tables.

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

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

1586 is True.

1587

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

1589 """

1590 name = 'Reference Band'

1591 shortname = 'refBand'

1592

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

1594

1595 @property

1596 def columns(self):

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

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

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

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

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

1602

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

1604 def getFilterAliasName(row):

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

1606 colName = row.idxmax()

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

1608

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

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

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

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

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

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

1615

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

1617 super().__init__(**kwargs)

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

1619

1620

1621class Photometry(Functor):

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

1623 # AB to NanoJansky (3631 Jansky).

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

1625 LOG_AB_FLUX_SCALE = 12.56

1626 FIVE_OVER_2LOG10 = 1.085736204758129569

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

1628 COADD_ZP = 27

1629

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

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

1632 self.col = colFlux

1633 self.colFluxErr = colFluxErr

1634

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

1636 self.fluxMag0Err = 0.

1637

1638 super().__init__(**kwargs)

1639

1640 @property

1641 def columns(self):

1642 return [self.col]

1643

1644 @property

1645 def name(self):

1646 return f'mag_{self.col}'

1647

1648 @classmethod

1649 def hypot(cls, a, b):

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

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

1652 a, b = b, a

1653 if a == 0.:

1654 return 0.

1655 q = b/a

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

1657

1658 def dn2flux(self, dn, fluxMag0):

1659 """Convert instrumental flux to nanojanskys."""

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

1661

1662 def dn2mag(self, dn, fluxMag0):

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

1664 with warnings.catch_warnings():

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

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

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

1668

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

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

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

1672 retVal *= self.AB_FLUX_SCALE / fluxMag0 / fluxMag0

1673 return retVal.astype(np.float32)

1674

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

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

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

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

1679

1680

1681class NanoJansky(Photometry):

1682 """Convert instrumental flux to nanojanskys."""

1683 def _func(self, df):

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

1685

1686

1687class NanoJanskyErr(Photometry):

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

1689 @property

1690 def columns(self):

1691 return [self.col, self.colFluxErr]

1692

1693 def _func(self, df):

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

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

1696

1697

1698class LocalPhotometry(Functor):

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

1700 the local photometric calibration.

1701

1702 Parameters

1703 ----------

1704 instFluxCol : `str`

1705 Name of the instrument flux column.

1706 instFluxErrCol : `str`

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

1708 photoCalibCol : `str`

1709 Name of local calibration column.

1710 photoCalibErrCol : `str`, optional

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

1712 be removed after v29.

1713

1714 See Also

1715 --------

1716 LocalNanojansky

1717 LocalNanojanskyErr

1718 """

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

1720

1721 def __init__(self,

1722 instFluxCol,

1723 instFluxErrCol,

1724 photoCalibCol,

1725 photoCalibErrCol=None,

1726 **kwargs):

1727 self.instFluxCol = instFluxCol

1728 self.instFluxErrCol = instFluxErrCol

1729 self.photoCalibCol = photoCalibCol

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

1731 if photoCalibErrCol is not None:

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

1733 category=FutureWarning)

1734 super().__init__(**kwargs)

1735

1736 def instFluxToNanojansky(self, instFlux, localCalib):

1737 """Convert instrument flux to nanojanskys.

1738

1739 Parameters

1740 ----------

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

1742 Array of instrument flux measurements.

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

1744 Array of local photometric calibration estimates.

1745

1746 Returns

1747 -------

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

1749 Array of calibrated flux measurements.

1750 """

1751 return instFlux * localCalib

1752

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

1754 """Convert instrument flux to nanojanskys.

1755

1756 Parameters

1757 ----------

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

1759 Array of instrument flux measurements. Ignored (accepted for

1760 backwards compatibility and consistency with magnitude-error

1761 calculation methods).

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

1763 Errors on associated ``instFlux`` values.

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

1765 Array of local photometric calibration estimates.

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

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

1768 will be removed after v29.

1769

1770 Returns

1771 -------

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

1773 Errors on calibrated flux measurements.

1774 """

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

1776 if localCalibErr is not None:

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

1778 category=FutureWarning)

1779 return instFluxErr * localCalib

1780

1781 def instFluxToMagnitude(self, instFlux, localCalib):

1782 """Convert instrument flux to nanojanskys.

1783

1784 Parameters

1785 ----------

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

1787 Array of instrument flux measurements.

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

1789 Array of local photometric calibration estimates.

1790

1791 Returns

1792 -------

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

1794 Array of calibrated AB magnitudes.

1795 """

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

1797

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

1799 """Convert instrument flux err to nanojanskys.

1800

1801 Parameters

1802 ----------

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

1804 Array of instrument flux measurements.

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

1806 Errors on associated ``instFlux`` values.

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

1808 Array of local photometric calibration estimates.

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

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

1811 will be removed after v29.

1812

1813 Returns

1814 -------

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

1816 Error on calibrated AB magnitudes.

1817 """

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

1819 if localCalibErr is not None:

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

1821 category=FutureWarning)

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

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

1824

1825

1826class LocalNanojansky(LocalPhotometry):

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

1828

1829 This returns units of nanojanskys.

1830 """

1831

1832 @property

1833 def columns(self):

1834 return [self.instFluxCol, self.photoCalibCol]

1835

1836 @property

1837 def name(self):

1838 return f'flux_{self.instFluxCol}'

1839

1840 def _func(self, df):

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

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

1843

1844

1845class LocalNanojanskyErr(LocalPhotometry):

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

1847

1848 This returns units of nanojanskys.

1849 """

1850

1851 @property

1852 def columns(self):

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

1854

1855 @property

1856 def name(self):

1857 return f'fluxErr_{self.instFluxCol}'

1858

1859 def _func(self, df):

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

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

1862

1863

1864class LocalDipoleMeanFlux(LocalPhotometry):

1865 """Compute absolute mean of dipole fluxes.

1866

1867 See Also

1868 --------

1869 LocalNanojansky

1870 LocalNanojanskyErr

1871 LocalDipoleMeanFluxErr

1872 LocalDipoleDiffFlux

1873 LocalDipoleDiffFluxErr

1874 """

1875 def __init__(self,

1876 instFluxPosCol,

1877 instFluxNegCol,

1878 instFluxPosErrCol,

1879 instFluxNegErrCol,

1880 photoCalibCol,

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

1882 photoCalibErrCol=None,

1883 **kwargs):

1884 self.instFluxNegCol = instFluxNegCol

1885 self.instFluxPosCol = instFluxPosCol

1886 self.instFluxNegErrCol = instFluxNegErrCol

1887 self.instFluxPosErrCol = instFluxPosErrCol

1888 self.photoCalibCol = photoCalibCol

1889 super().__init__(instFluxNegCol,

1890 instFluxNegErrCol,

1891 photoCalibCol,

1892 photoCalibErrCol,

1893 **kwargs)

1894

1895 @property

1896 def columns(self):

1897 return [self.instFluxPosCol,

1898 self.instFluxNegCol,

1899 self.photoCalibCol]

1900

1901 @property

1902 def name(self):

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

1904

1905 def _func(self, df):

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

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

1908

1909

1910class LocalDipoleMeanFluxErr(LocalDipoleMeanFlux):

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

1912

1913 See Also

1914 --------

1915 LocalNanojansky

1916 LocalNanojanskyErr

1917 LocalDipoleMeanFlux

1918 LocalDipoleDiffFlux

1919 LocalDipoleDiffFluxErr

1920 """

1921

1922 @property

1923 def columns(self):

1924 return [self.instFluxPosCol,

1925 self.instFluxNegCol,

1926 self.instFluxPosErrCol,

1927 self.instFluxNegErrCol,

1928 self.photoCalibCol]

1929

1930 @property

1931 def name(self):

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

1933

1934 def _func(self, df):

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

1936

1937

1938class LocalDipoleDiffFlux(LocalDipoleMeanFlux):

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

1940

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

1942

1943 See Also

1944 --------

1945 LocalNanojansky

1946 LocalNanojanskyErr

1947 LocalDipoleMeanFlux

1948 LocalDipoleMeanFluxErr

1949 LocalDipoleDiffFluxErr

1950 """

1951

1952 @property

1953 def columns(self):

1954 return [self.instFluxPosCol,

1955 self.instFluxNegCol,

1956 self.photoCalibCol]

1957

1958 @property

1959 def name(self):

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

1961

1962 def _func(self, df):

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

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

1965

1966

1967class LocalDipoleDiffFluxErr(LocalDipoleMeanFlux):

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

1969

1970 See Also

1971 --------

1972 LocalNanojansky

1973 LocalNanojanskyErr

1974 LocalDipoleMeanFlux

1975 LocalDipoleMeanFluxErr

1976 LocalDipoleDiffFlux

1977 """

1978

1979 @property

1980 def columns(self):

1981 return [self.instFluxPosCol,

1982 self.instFluxNegCol,

1983 self.instFluxPosErrCol,

1984 self.instFluxNegErrCol,

1985 self.photoCalibCol]

1986

1987 @property

1988 def name(self):

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

1990

1991 def _func(self, df):

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

1993

1994

1995class Ebv(Functor):

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

1997 _defaultDataset = 'ref'

1998 name = "E(B-V)"

1999 shortname = "ebv"

2000

2001 def __init__(self, **kwargs):

2002 # Import is only needed for Ebv.

2003 # Suppress unnecessary .dustmapsrc log message on import.

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

2005 with redirect_stdout(devnull):

2006 from dustmaps.sfd import SFDQuery

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

2008 self.sfd = SFDQuery()

2009 super().__init__(**kwargs)

2010

2011 def _func(self, df):

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

2013 ebv = self.sfd(coords)

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

2015

2016

2017class MomentsBase(Functor):

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

2019

2020 Attributes

2021 ----------

2022 is_covariance : bool

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

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

2025 """

2026

2027 is_covariance: bool = True

2028

2029 def __init__(self,

2030 shape_1_1,

2031 shape_2_2,

2032 shape_1_2,

2033 colCD_1_1,

2034 colCD_1_2,

2035 colCD_2_1,

2036 colCD_2_2,

2037 **kwargs):

2038 self.shape_1_1 = shape_1_1

2039 self.shape_2_2 = shape_2_2

2040 self.shape_1_2 = shape_1_2

2041 self.colCD_1_1 = colCD_1_1

2042 self.colCD_1_2 = colCD_1_2

2043 self.colCD_2_1 = colCD_2_1

2044 self.colCD_2_2 = colCD_2_2

2045 super().__init__(**kwargs)

2046

2047 @property

2048 def columns(self):

2049 return [

2050 self.shape_1_1,

2051 self.shape_2_2,

2052 self.shape_1_2,

2053 ] + self.columns_ref

2054

2055 @property

2056 def columns_ref(self):

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

2058 return [

2059 self.colCD_1_1,

2060 self.colCD_1_2,

2061 self.colCD_2_1,

2062 self.colCD_2_2]

2063

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

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

2066

2067 Parameters

2068 ----------

2069 df

2070 The data frame.

2071 sky

2072 Whether to compute the terms in sky coordinates.

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

2074 UU, VV and UV.

2075

2076 Returns

2077 -------

2078 xx_p_yy

2079 The sum of the diagonal terms of the covariance.

2080 xx_m_yy

2081 The difference of the diagonal terms of the covariance.

2082 t2

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

2084 """

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

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

2087 xx_m_yy = xx - yy

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

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

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

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

2092 return xx + yy, xx_m_yy, t2

2093

2094 def get_xx(self, df):

2095 xx = df[self.shape_1_1]

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

2097

2098 def get_yy(self, df):

2099 yy = df[self.shape_2_2]

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

2101

2102 def get_xy(self, df):

2103 xy = df[self.shape_1_2]

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

2105

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

2107 # CD_matrix * moments_matrix * CD_matrix.T

2108 def sky_uu(self, df):

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

2110 i_xx = self.get_xx(df)

2111 i_yy = self.get_yy(df)

2112 i_xy = self.get_xy(df)

2113 CD_1_1 = df[self.colCD_1_1]

2114 CD_1_2 = df[self.colCD_1_2]

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

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

2117

2118 def sky_vv(self, df):

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

2120 i_xx = self.get_xx(df)

2121 i_yy = self.get_yy(df)

2122 i_xy = self.get_xy(df)

2123 CD_2_1 = df[self.colCD_2_1]

2124 CD_2_2 = df[self.colCD_2_2]

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

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

2127

2128 def sky_uv(self, df):

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

2130 i_xx = self.get_xx(df)

2131 i_yy = self.get_yy(df)

2132 i_xy = self.get_xy(df)

2133 CD_1_1 = df[self.colCD_1_1]

2134 CD_1_2 = df[self.colCD_1_2]

2135 CD_2_1 = df[self.colCD_2_1]

2136 CD_2_2 = df[self.colCD_2_2]

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

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

2139

2140 def get_g1(self, df):

2141 """

2142 Calculate shear-type ellipticity parameter G1.

2143 """

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

2145 sky_uu = self.sky_uu(df)

2146 sky_vv = self.sky_vv(df)

2147 sky_uv = self.sky_uv(df)

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

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

2150

2151 def get_g2(self, df):

2152 """

2153 Calculate shear-type ellipticity parameter G2.

2154

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

2156 sign convention.

2157 """

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

2159 sky_uu = self.sky_uu(df)

2160 sky_vv = self.sky_vv(df)

2161 sky_uv = self.sky_uv(df)

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

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

2164

2165 def get_trace(self, df):

2166 sky_uu = self.sky_uu(df)

2167 sky_vv = self.sky_vv(df)

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

2169

2170

2171class MomentsG1Sky(MomentsBase):

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

2173 shear parameterization"""

2174 _defaultDataset = 'meas'

2175 name = "moments_g1"

2176 shortname = "moments_g1"

2177

2178 def _func(self, df):

2179 sky_g1 = self.get_g1(df)

2180

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

2182

2183

2184class MomentsG2Sky(MomentsBase):

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

2186 shear parameterization"""

2187 _defaultDataset = 'meas'

2188 name = "moments_g2"

2189 shortname = "moments_g2"

2190

2191 def _func(self, df):

2192 sky_g2 = self.get_g2(df)

2193

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

2195

2196

2197class MomentsTraceSky(MomentsBase):

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

2199

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

2201 half of the trace of the second moments tensor.

2202 """

2203 _defaultDataset = 'meas'

2204 name = "moments_trace"

2205 shortname = "moments_trace"

2206

2207 def _func(self, df):

2208 sky_trace_radians = self.get_trace(df)

2209

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

2211

2212

2213class MomentsIuuSky(MomentsBase):

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

2215 _defaultDataset = 'meas'

2216 name = "moments_uu"

2217 shortname = "moments_uu"

2218

2219 def _func(self, df):

2220 sky_uu_radians = self.sky_uu(df)

2221

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

2223

2224

2225class CorrelationIuuSky(MomentsIuuSky):

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

2227 is_covariance = False

2228

2229

2230class MomentsIvvSky(MomentsBase):

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

2232 _defaultDataset = 'meas'

2233 name = "moments_vv"

2234 shortname = "moments_vv"

2235

2236 def _func(self, df):

2237 sky_vv_radians = self.sky_vv(df)

2238

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

2240

2241

2242class CorrelationIvvSky(MomentsIvvSky):

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

2244 is_covariance = False

2245

2246

2247class MomentsIuvSky(MomentsBase):

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

2249 _defaultDataset = 'meas'

2250 name = "moments_uv"

2251 shortname = "moments_uv"

2252

2253 def _func(self, df):

2254 sky_uv_radians = self.sky_uv(df)

2255

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

2257

2258

2259class CorrelationIuvSky(MomentsIuvSky):

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

2261 is_covariance = False

2262

2263

2264class PositionAngleFromMoments(MomentsBase):

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

2266 _defaultDataset = 'meas'

2267 name = "moments_theta"

2268 shortname = "moments_theta"

2269

2270 def _func(self, df):

2271 sky_uu = self.sky_uu(df)

2272 sky_vv = self.sky_vv(df)

2273 sky_uv = self.sky_uv(df)

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

2275

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

2277

2278

2279class PositionAngleFromCorrelation(PositionAngleFromMoments):

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

2281 is_covariance = False

2282

2283

2284class SemimajorAxisFromMoments(MomentsBase):

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

2286 _defaultDataset = 'meas'

2287 name = "moments_a"

2288 shortname = "moments_a"

2289

2290 def _func(self, df):

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

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

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

2294

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

2296

2297

2298class SemimajorAxisFromCorrelation(SemimajorAxisFromMoments):

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

2300 is_covariance = False

2301

2302

2303class SemiminorAxisFromMoments(MomentsBase):

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

2305 _defaultDataset = 'meas'

2306 name = "moments_b"

2307 shortname = "moments_b"

2308

2309 def _func(self, df):

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

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

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

2313

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

2315

2316

2317class SemiminorAxisFromCorrelation(SemiminorAxisFromMoments):

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

2319 is_covariance = False