Coverage for python/lsst/pipe/tasks/functors.py: 86%
951 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 09:36 +0000
1# This file is part of pipe_tasks.
2#
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.
8#
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.resources import ResourcePath
62from lsst.utils import doImport
63from lsst.utils.introspection import get_full_type_name
66def init_fromDict(initDict, basePath='lsst.pipe.tasks.functors',
67 typeKey='functor', name=None):
68 """Initialize an object defined in a dictionary.
70 The object needs to be importable as f'{basePath}.{initDict[typeKey]}'.
71 The positional and keyword arguments (if any) are contained in "args" and
72 "kwargs" entries in the dictionary, respectively.
73 This is used in `~lsst.pipe.tasks.functors.CompositeFunctor.from_yaml` to
74 initialize a composite functor from a specification in a YAML file.
76 Parameters
77 ----------
78 initDict : dictionary
79 Dictionary describing object's initialization.
80 Must contain an entry keyed by ``typeKey`` that is the name of the
81 object, relative to ``basePath``.
82 basePath : str
83 Path relative to module in which ``initDict[typeKey]`` is defined.
84 typeKey : str
85 Key of ``initDict`` that is the name of the object (relative to
86 ``basePath``).
87 """
88 initDict = initDict.copy()
89 # TO DO: DM-21956 We should be able to define functors outside this module
90 pythonType = doImport(f'{basePath}.{initDict.pop(typeKey)}')
91 args = []
92 if 'args' in initDict:
93 args = initDict.pop('args')
94 if isinstance(args, str):
95 args = [args]
96 try:
97 element = pythonType(*args, **initDict)
98 except Exception as e:
99 message = f'Error in constructing functor "{name}" of type {pythonType.__name__} with args: {args}'
100 raise type(e)(message, e.args)
101 return element
104class Functor(object):
105 """Define and execute a calculation on a DataFrame or Handle holding a
106 DataFrame.
108 The `__call__` method accepts either a `~pandas.DataFrame` object or a
109 `~lsst.daf.butler.DeferredDatasetHandle` or
110 `~lsst.pipe.base.InMemoryDatasetHandle`, and returns the
111 result of the calculation as a single column.
112 Each functor defines what columns are needed for the calculation, and only
113 these columns are read from the dataset handle.
115 The action of `__call__` consists of two steps: first, loading the
116 necessary columns from disk into memory as a `~pandas.DataFrame` object;
117 and second, performing the computation on this DataFrame and returning the
118 result.
120 To define a new `Functor`, a subclass must define a `_func` method,
121 that takes a `~pandas.DataFrame` and returns result in a `~pandas.Series`.
122 In addition, it must define the following attributes:
124 * `_columns`: The columns necessary to perform the calculation
125 * `name`: A name appropriate for a figure axis label
126 * `shortname`: A name appropriate for use as a dictionary key
128 On initialization, a `Functor` should declare what band (``filt`` kwarg)
129 and dataset (e.g. ``'ref'``, ``'meas'``, ``'forced_src'``) it is intended
130 to be applied to.
131 This enables the `_get_data` method to extract the proper columns from the
132 underlying data.
133 If not specified, the dataset will fall back on the `_defaultDataset`
134 attribute.
135 If band is not specified and ``dataset`` is anything other than ``'ref'``,
136 then an error will be raised when trying to perform the calculation.
138 Originally, `Functor` was set up to expect datasets formatted like the
139 ``deepCoadd_obj`` dataset; that is, a DataFrame with a multi-level column
140 index, with the levels of the column index being ``band``, ``dataset``, and
141 ``column``.
142 It has since been generalized to apply to DataFrames without multi-level
143 indices and multi-level indices with just ``dataset`` and ``column``
144 levels.
145 In addition, the `_get_data` method that reads the columns from the
146 underlying data will return a DataFrame with column index levels defined by
147 the `_dfLevels` attribute; by default, this is ``column``.
149 The `_dfLevels` attributes should generally not need to be changed, unless
150 `_func` needs columns from multiple filters or datasets to do the
151 calculation.
152 An example of this is the `~lsst.pipe.tasks.functors.Color` functor, for
153 which `_dfLevels = ('band', 'column')`, and `_func` expects the DataFrame
154 it gets to have those levels in the column index.
156 Parameters
157 ----------
158 filt : str
159 Band upon which to do the calculation.
161 dataset : str
162 Dataset upon which to do the calculation (e.g., 'ref', 'meas',
163 'forced_src').
164 """
166 _defaultDataset = 'ref'
167 _dfLevels = ('column',)
168 _defaultNoDup = False
170 def __init__(self, filt=None, dataset=None, noDup=None):
171 self.filt = filt
172 self.dataset = dataset if dataset is not None else self._defaultDataset
173 self._noDup = noDup
174 self.log = logging.getLogger(type(self).__name__)
176 @property
177 def noDup(self):
178 """Do not explode by band if used on object table."""
179 if self._noDup is not None: 179 ↛ 180line 179 didn't jump to line 180 because the condition on line 179 was never true
180 return self._noDup
181 else:
182 return self._defaultNoDup
184 @property
185 def columns(self):
186 """Columns required to perform calculation."""
187 if not hasattr(self, '_columns'): 187 ↛ 188line 187 didn't jump to line 188 because the condition on line 187 was never true
188 raise NotImplementedError('Must define columns property or _columns attribute')
189 return self._columns
191 def _get_data_columnLevels(self, data, columnIndex=None):
192 """Gets the names of the column index levels.
194 This should only be called in the context of a multilevel table.
196 Parameters
197 ----------
198 data : various
199 The data to be read, can be a
200 `~lsst.daf.butler.DeferredDatasetHandle` or
201 `~lsst.pipe.base.InMemoryDatasetHandle`.
202 columnIndex (optional): pandas `~pandas.Index` object
203 If not passed, then it is read from the
204 `~lsst.daf.butler.DeferredDatasetHandle`
205 for `~lsst.pipe.base.InMemoryDatasetHandle`.
206 """
207 if columnIndex is None: 207 ↛ 208line 207 didn't jump to line 208 because the condition on line 207 was never true
208 columnIndex = data.get(component="columns")
209 return columnIndex.names
211 def _get_data_columnLevelNames(self, data, columnIndex=None):
212 """Gets the content of each of the column levels for a multilevel
213 table.
214 """
215 if columnIndex is None: 215 ↛ 216line 215 didn't jump to line 216 because the condition on line 215 was never true
216 columnIndex = data.get(component="columns")
218 columnLevels = columnIndex.names
219 columnLevelNames = {
220 level: list(np.unique(np.array([c for c in columnIndex])[:, i]))
221 for i, level in enumerate(columnLevels)
222 }
223 return columnLevelNames
225 def _colsFromDict(self, colDict, columnIndex=None):
226 """Converts dictionary column specficiation to a list of columns."""
227 new_colDict = {}
228 columnLevels = self._get_data_columnLevels(None, columnIndex=columnIndex)
230 for i, lev in enumerate(columnLevels):
231 if lev in colDict: 231 ↛ 237line 231 didn't jump to line 237 because the condition on line 231 was always true
232 if isinstance(colDict[lev], str):
233 new_colDict[lev] = [colDict[lev]]
234 else:
235 new_colDict[lev] = colDict[lev]
236 else:
237 new_colDict[lev] = columnIndex.levels[i]
239 levelCols = [new_colDict[lev] for lev in columnLevels]
240 cols = list(product(*levelCols))
241 colsAvailable = [col for col in cols if col in columnIndex]
242 return colsAvailable
244 def multilevelColumns(self, data, columnIndex=None, returnTuple=False):
245 """Returns columns needed by functor from multilevel dataset.
247 To access tables with multilevel column structure, the
248 `~lsst.daf.butler.DeferredDatasetHandle` or
249 `~lsst.pipe.base.InMemoryDatasetHandle` needs to be passed
250 either a list of tuples or a dictionary.
252 Parameters
253 ----------
254 data : various
255 The data as either `~lsst.daf.butler.DeferredDatasetHandle`, or
256 `~lsst.pipe.base.InMemoryDatasetHandle`.
257 columnIndex (optional): pandas `~pandas.Index` object
258 Either passed or read in from
259 `~lsst.daf.butler.DeferredDatasetHandle`.
260 `returnTuple` : `bool`
261 If true, then return a list of tuples rather than the column
262 dictionary specification.
263 This is set to `True` by `CompositeFunctor` in order to be able to
264 combine columns from the various component functors.
266 """
267 if not isinstance(data, (DeferredDatasetHandle, InMemoryDatasetHandle)): 267 ↛ 268line 267 didn't jump to line 268 because the condition on line 267 was never true
268 raise RuntimeError(f"Unexpected data type. Got {get_full_type_name(data)}.")
270 if columnIndex is None: 270 ↛ 271line 270 didn't jump to line 271 because the condition on line 270 was never true
271 columnIndex = data.get(component="columns")
273 # Confirm that the dataset has the column levels the functor is
274 # expecting it to have.
275 columnLevels = self._get_data_columnLevels(data, columnIndex)
277 columnDict = {'column': self.columns,
278 'dataset': self.dataset}
279 if self.filt is None:
280 columnLevelNames = self._get_data_columnLevelNames(data, columnIndex)
281 if "band" in columnLevels: 281 ↛ 293line 281 didn't jump to line 293 because the condition on line 281 was always true
282 if self.dataset == "ref": 282 ↛ 285line 282 didn't jump to line 285 because the condition on line 282 was always true
283 columnDict["band"] = columnLevelNames["band"][0]
284 else:
285 raise ValueError(f"'filt' not set for functor {self.name}"
286 f"(dataset {self.dataset}) "
287 "and DataFrame "
288 "contains multiple filters in column index. "
289 "Set 'filt' or set 'dataset' to 'ref'.")
290 else:
291 columnDict['band'] = self.filt
293 if returnTuple:
294 return self._colsFromDict(columnDict, columnIndex=columnIndex)
295 else:
296 return columnDict
298 def _func(self, df, dropna=True):
299 raise NotImplementedError('Must define calculation on DataFrame')
301 def _get_columnIndex(self, data):
302 """Return columnIndex."""
304 if isinstance(data, (DeferredDatasetHandle, InMemoryDatasetHandle)): 304 ↛ 307line 304 didn't jump to line 307 because the condition on line 304 was always true
305 return data.get(component="columns")
306 else:
307 return None
309 def _get_data(self, data):
310 """Retrieve DataFrame necessary for calculation.
312 The data argument can be a `~pandas.DataFrame`, a
313 `~lsst.daf.butler.DeferredDatasetHandle`, or
314 an `~lsst.pipe.base.InMemoryDatasetHandle`.
316 Returns a DataFrame upon which `self._func` can act.
317 """
318 # We wrap a DataFrame in a handle here to take advantage of the
319 # DataFrame delegate DataFrame column wrangling abilities.
320 if isinstance(data, pd.DataFrame):
321 _data = InMemoryDatasetHandle(data, storageClass="DataFrame")
322 elif isinstance(data, (DeferredDatasetHandle, InMemoryDatasetHandle)): 322 ↛ 325line 322 didn't jump to line 325 because the condition on line 322 was always true
323 _data = data
324 else:
325 raise RuntimeError(f"Unexpected type provided for data. Got {get_full_type_name(data)}.")
327 # First thing to do: check to see if the data source has a multilevel
328 # column index or not.
329 columnIndex = self._get_columnIndex(_data)
330 is_multiLevel = isinstance(columnIndex, pd.MultiIndex)
332 # Get proper columns specification for this functor.
333 if is_multiLevel:
334 columns = self.multilevelColumns(_data, columnIndex=columnIndex)
335 else:
336 columns = self.columns
338 # Load in-memory DataFrame with appropriate columns the gen3 way.
339 df = _data.get(parameters={"columns": columns})
341 # Drop unnecessary column levels.
342 if is_multiLevel:
343 df = self._setLevels(df)
345 return df
347 def _setLevels(self, df):
348 levelsToDrop = [n for n in df.columns.names if n not in self._dfLevels]
349 df.columns = df.columns.droplevel(levelsToDrop)
350 return df
352 def _dropna(self, vals):
353 return vals.dropna()
355 def __call__(self, data, dropna=False):
356 df = self._get_data(data)
357 try:
358 vals = self._func(df)
359 except Exception as e:
360 self.log.error("Exception in %s call: %s: %s", self.name, type(e).__name__, e)
361 vals = self.fail(df)
362 if dropna:
363 vals = self._dropna(vals)
365 return vals
367 def difference(self, data1, data2, **kwargs):
368 """Computes difference between functor called on two different
369 DataFrame/Handle objects.
370 """
371 return self(data1, **kwargs) - self(data2, **kwargs)
373 def fail(self, df):
374 return pd.Series(np.full(len(df), np.nan), index=df.index)
376 @property
377 def name(self):
378 """Full name of functor (suitable for figure labels)."""
379 return NotImplementedError
381 @property
382 def shortname(self):
383 """Short name of functor (suitable for column name/dict key)."""
384 return self.name
387class CompositeFunctor(Functor):
388 """Perform multiple calculations at once on a catalog.
390 The role of a `CompositeFunctor` is to group together computations from
391 multiple functors.
392 Instead of returning `~pandas.Series` a `CompositeFunctor` returns a
393 `~pandas.DataFrame`, with the column names being the keys of ``funcDict``.
395 The `columns` attribute of a `CompositeFunctor` is the union of all columns
396 in all the component functors.
398 A `CompositeFunctor` does not use a `_func` method itself; rather, when a
399 `CompositeFunctor` is called, all its columns are loaded at once, and the
400 resulting DataFrame is passed to the `_func` method of each component
401 functor.
402 This has the advantage of only doing I/O (reading from parquet file) once,
403 and works because each individual `_func` method of each component functor
404 does not care if there are *extra* columns in the DataFrame being passed;
405 only that it must contain *at least* the `columns` it expects.
407 An important and useful class method is `from_yaml`, which takes as an
408 argument the path to a YAML file specifying a collection of functors.
410 Parameters
411 ----------
412 funcs : `dict` or `list`
413 Dictionary or list of functors.
414 If a list, then it will be converted into a dictonary according to the
415 `.shortname` attribute of each functor.
416 """
417 dataset = None
418 name = "CompositeFunctor"
420 def __init__(self, funcs, **kwargs):
422 if type(funcs) is dict:
423 self.funcDict = funcs
424 else:
425 self.funcDict = {f.shortname: f for f in funcs}
427 self._filt = None
429 super().__init__(**kwargs)
431 @property
432 def filt(self):
433 return self._filt
435 @filt.setter
436 def filt(self, filt):
437 if filt is not None:
438 for _, f in self.funcDict.items():
439 f.filt = filt
440 self._filt = filt
442 def update(self, new):
443 """Update the functor with new functors."""
444 if isinstance(new, dict): 444 ↛ 446line 444 didn't jump to line 446 because the condition on line 444 was always true
445 self.funcDict.update(new)
446 elif isinstance(new, CompositeFunctor):
447 self.funcDict.update(new.funcDict)
448 else:
449 raise TypeError('Can only update with dictionary or CompositeFunctor.')
451 # Make sure new functors have the same 'filt' set.
452 if self.filt is not None: 452 ↛ 453line 452 didn't jump to line 453 because the condition on line 452 was never true
453 self.filt = self.filt
455 @property
456 def columns(self):
457 return list(set([x for y in [f.columns for f in self.funcDict.values()] for x in y]))
459 def multilevelColumns(self, data, **kwargs):
460 # Get the union of columns for all component functors.
461 # Note the need to have `returnTuple=True` here.
462 return list(
463 set(
464 [
465 x
466 for y in [
467 f.multilevelColumns(data, returnTuple=True, **kwargs) for f in self.funcDict.values()
468 ]
469 for x in y
470 ]
471 )
472 )
474 def __call__(self, data, **kwargs):
475 """Apply the functor to the data table.
477 Parameters
478 ----------
479 data : various
480 The data represented as `~lsst.daf.butler.DeferredDatasetHandle`,
481 `~lsst.pipe.base.InMemoryDatasetHandle`, or `~pandas.DataFrame`.
482 The table or a pointer to a table on disk from which columns can
483 be accessed.
484 """
485 if isinstance(data, pd.DataFrame):
486 _data = InMemoryDatasetHandle(data, storageClass="DataFrame")
487 elif isinstance(data, (DeferredDatasetHandle, InMemoryDatasetHandle)): 487 ↛ 490line 487 didn't jump to line 490 because the condition on line 487 was always true
488 _data = data
489 else:
490 raise RuntimeError(f"Unexpected type provided for data. Got {get_full_type_name(data)}.")
492 columnIndex = self._get_columnIndex(_data)
494 if isinstance(columnIndex, pd.MultiIndex):
495 columns = self.multilevelColumns(_data, columnIndex=columnIndex)
496 df = _data.get(parameters={"columns": columns})
498 valDict = {}
499 for k, f in self.funcDict.items():
500 try:
501 subdf = f._setLevels(
502 df[f.multilevelColumns(_data, returnTuple=True, columnIndex=columnIndex)]
503 )
504 valDict[k] = f._func(subdf)
505 except Exception as e:
506 self.log.exception(
507 "Exception in %s (funcs: %s) call: %s",
508 self.name,
509 str(list(self.funcDict.keys())),
510 type(e).__name__,
511 )
512 try:
513 valDict[k] = f.fail(subdf)
514 except NameError:
515 raise e
517 else:
518 df = _data.get(parameters={"columns": self.columns})
520 valDict = {k: f._func(df) for k, f in self.funcDict.items()}
522 # Check that output columns are actually columns.
523 for name, colVal in valDict.items():
524 if len(colVal.shape) != 1: 524 ↛ 525line 524 didn't jump to line 525 because the condition on line 524 was never true
525 raise RuntimeError("Transformed column '%s' is not the shape of a column. "
526 "It is shaped %s and type %s." % (name, colVal.shape, type(colVal)))
528 try:
529 valDf = pd.concat(valDict, axis=1)
530 except TypeError:
531 print([(k, type(v)) for k, v in valDict.items()])
532 raise
534 if kwargs.get('dropna', False):
535 valDf = valDf.dropna(how='any')
537 return valDf
539 @classmethod
540 def renameCol(cls, col, renameRules):
541 if renameRules is None: 541 ↛ 542line 541 didn't jump to line 542 because the condition on line 541 was never true
542 return col
543 for old, new in renameRules:
544 if col.startswith(old):
545 col = col.replace(old, new)
546 return col
548 @classmethod
549 def from_file(cls, filename, **kwargs):
550 # Allow environment variables in the filename.
551 filename = os.path.expandvars(filename)
552 with ResourcePath(filename).open("r") as f:
553 translationDefinition = yaml.safe_load(f)
555 return cls.from_yaml(translationDefinition, **kwargs)
557 @classmethod
558 def from_yaml(cls, translationDefinition, **kwargs):
559 funcs = {}
560 for func, val in translationDefinition['funcs'].items():
561 funcs[func] = init_fromDict(val, name=func)
563 if 'flag_rename_rules' in translationDefinition: 563 ↛ 566line 563 didn't jump to line 566 because the condition on line 563 was always true
564 renameRules = translationDefinition['flag_rename_rules']
565 else:
566 renameRules = None
568 if 'calexpFlags' in translationDefinition: 568 ↛ 569line 568 didn't jump to line 569 because the condition on line 568 was never true
569 for flag in translationDefinition['calexpFlags']:
570 funcs[cls.renameCol(flag, renameRules)] = Column(flag, dataset='calexp')
572 if 'refFlags' in translationDefinition:
573 for flag in translationDefinition['refFlags']:
574 funcs[cls.renameCol(flag, renameRules)] = Column(flag, dataset='ref')
576 if 'forcedFlags' in translationDefinition:
577 for flag in translationDefinition['forcedFlags']:
578 funcs[cls.renameCol(flag, renameRules)] = Column(flag, dataset='forced_src')
580 if 'flags' in translationDefinition: 580 ↛ 584line 580 didn't jump to line 584 because the condition on line 580 was always true
581 for flag in translationDefinition['flags']:
582 funcs[cls.renameCol(flag, renameRules)] = Column(flag, dataset='meas')
584 return cls(funcs, **kwargs)
587def mag_aware_eval(df, expr, log):
588 """Evaluate an expression on a DataFrame, knowing what the 'mag' function
589 means.
591 Builds on `pandas.DataFrame.eval`, which parses and executes math on
592 DataFrames.
594 Parameters
595 ----------
596 df : ~pandas.DataFrame
597 DataFrame on which to evaluate expression.
599 expr : str
600 Expression.
601 """
602 try:
603 expr_new = re.sub(r'mag\((\w+)\)', r'-2.5*log(\g<1>)/log(10)', expr)
604 val = df.eval(expr_new)
605 except Exception as e: # Should check what actually gets raised
606 log.error("Exception in mag_aware_eval: %s: %s", type(e).__name__, e)
607 expr_new = re.sub(r'mag\((\w+)\)', r'-2.5*log(\g<1>_instFlux)/log(10)', expr)
608 val = df.eval(expr_new)
609 return val
612class CustomFunctor(Functor):
613 """Arbitrary computation on a catalog.
615 Column names (and thus the columns to be loaded from catalog) are found by
616 finding all words and trying to ignore all "math-y" words.
618 Parameters
619 ----------
620 expr : str
621 Expression to evaluate, to be parsed and executed by
622 `~lsst.pipe.tasks.functors.mag_aware_eval`.
623 """
624 _ignore_words = ('mag', 'sin', 'cos', 'exp', 'log', 'sqrt')
626 def __init__(self, expr, **kwargs):
627 self.expr = expr
628 super().__init__(**kwargs)
630 @property
631 def name(self):
632 return self.expr
634 @property
635 def columns(self):
636 flux_cols = re.findall(r'mag\(\s*(\w+)\s*\)', self.expr)
638 cols = [c for c in re.findall(r'[a-zA-Z_]+', self.expr) if c not in self._ignore_words]
639 not_a_col = []
640 for c in flux_cols: 640 ↛ 641line 640 didn't jump to line 641 because the loop on line 640 never started
641 if not re.search('_instFlux$', c):
642 cols.append(f'{c}_instFlux')
643 not_a_col.append(c)
644 else:
645 cols.append(c)
647 return list(set([c for c in cols if c not in not_a_col]))
649 def _func(self, df):
650 return mag_aware_eval(df, self.expr, self.log)
653class Column(Functor):
654 """Get column with a specified name."""
656 def __init__(self, col, **kwargs):
657 self.col = col
658 super().__init__(**kwargs)
660 @property
661 def name(self):
662 return self.col
664 @property
665 def columns(self):
666 return [self.col]
668 def _func(self, df):
669 return df[self.col]
672class Index(Functor):
673 """Return the value of the index for each object."""
675 columns = ['coord_ra'] # Just a dummy; something has to be here.
676 _defaultDataset = 'ref'
677 _defaultNoDup = True
679 def _func(self, df):
680 return pd.Series(df.index, index=df.index)
683class CoordColumn(Column):
684 """Base class for coordinate column, in degrees."""
685 _radians = True
687 def __init__(self, col, **kwargs):
688 super().__init__(col, **kwargs)
690 def _func(self, df):
691 # Must not modify original column in case that column is used by
692 # another functor.
693 output = df[self.col] * 180 / np.pi if self._radians else df[self.col]
694 return output
697class RAColumn(CoordColumn):
698 """Right Ascension, in degrees."""
699 name = 'RA'
700 _defaultNoDup = True
702 def __init__(self, **kwargs):
703 super().__init__('coord_ra', **kwargs)
705 def __call__(self, catalog, **kwargs):
706 return super().__call__(catalog, **kwargs)
709class DecColumn(CoordColumn):
710 """Declination, in degrees."""
711 name = 'Dec'
712 _defaultNoDup = True
714 def __init__(self, **kwargs):
715 super().__init__('coord_dec', **kwargs)
717 def __call__(self, catalog, **kwargs):
718 return super().__call__(catalog, **kwargs)
721class RAErrColumn(CoordColumn):
722 """Angular uncertainty in Right Ascension, in degrees.
724 Note that ``coord_raErr`` is the uncertainty on RA·cos(Dec),
725 i.e. an arc-length on the sky in the RA direction.
726 """
727 name = 'RAErr'
728 _defaultNoDup = True
730 def __init__(self, **kwargs):
731 super().__init__('coord_raErr', **kwargs)
734class DecErrColumn(CoordColumn):
735 """Uncertainty in declination, in degrees."""
736 name = 'DecErr'
737 _defaultNoDup = True
739 def __init__(self, **kwargs):
740 super().__init__('coord_decErr', **kwargs)
743class RADecCovColumn(Column):
744 """Tangent-plane angular covariance, in degrees^2.
746 As with the RA error, this is the covariance between RA·cos(Dec) and Dec,
747 Cov(xi, eta) = Cov(RA*cos(Dec), Dec).
748 """
749 _radians = True
750 name = 'RADecCov'
751 _defaultNoDup = True
753 def __init__(self, **kwargs):
754 super().__init__('coord_ra_dec_Cov', **kwargs)
756 def _func(self, df):
757 # Must not modify original column in case that column is used by
758 # another functor.
759 output = df[self.col]*(180/np.pi)**2 if self._radians else df[self.col]
760 return output
763class MultibandColumn(Column):
764 """A column with a band in a multiband table."""
765 def __init__(self, col, band_to_check, **kwargs):
766 self._band_to_check = band_to_check
767 super().__init__(col=col, **kwargs)
769 @property
770 def band_to_check(self):
771 return self._band_to_check
774class MultibandSinglePrecisionFloatColumn(MultibandColumn):
775 """A float32 MultibandColumn"""
776 def _func(self, df):
777 return super()._func(df).astype(np.float32)
780class SinglePrecisionFloatColumn(Column):
781 """Return a column cast to a single-precision float."""
783 def _func(self, df):
784 return df[self.col].astype(np.float32)
787class HtmIndex20(Functor):
788 """Compute the level 20 HtmIndex for the catalog.
790 Notes
791 -----
792 This functor was implemented to satisfy requirements of old APDB interface
793 which required the ``pixelId`` column in DiaObject with HTM20 index.
794 The APDB interface had migrated to not need that information, but we keep
795 this class in case it may be useful for something else.
796 """
797 name = "Htm20"
798 htmLevel = 20
799 _radians = True
801 def __init__(self, ra, dec, **kwargs):
802 self.pixelator = sphgeom.HtmPixelization(self.htmLevel)
803 self.ra = ra
804 self.dec = dec
805 self._columns = [self.ra, self.dec]
806 super().__init__(**kwargs)
808 def _func(self, df):
810 def computePixel(row):
811 if self._radians: 811 ↛ 816line 811 didn't jump to line 816 because the condition on line 811 was always true
812 sphPoint = geom.SpherePoint(row[self.ra],
813 row[self.dec],
814 geom.radians)
815 else:
816 sphPoint = geom.SpherePoint(row[self.ra],
817 row[self.dec],
818 geom.degrees)
819 return self.pixelator.index(sphPoint.getVector())
821 return df.apply(computePixel, axis=1, result_type='reduce').astype('int64')
824def fluxName(col):
825 """Append _instFlux to the column name if it doesn't have it already."""
826 if not col.endswith('_instFlux'): 826 ↛ 828line 826 didn't jump to line 828 because the condition on line 826 was always true
827 col += '_instFlux'
828 return col
831def fluxErrName(col):
832 """Append _instFluxErr to the column name if it doesn't have it already."""
833 if not col.endswith('_instFluxErr'):
834 col += '_instFluxErr'
835 return col
838class Mag(Functor):
839 """Compute calibrated magnitude.
841 Returns the flux at mag=0.
842 The default ``fluxMag0`` is 63095734448.0194, which is default for HSC.
843 TO DO: This default should be made configurable in DM-21955.
845 This calculation hides warnings about invalid values and dividing by zero.
847 As with all functors, a ``dataset`` and ``filt`` kwarg should be provided
848 upon initialization.
849 Unlike the default `Functor`, however, the default dataset for a `Mag` is
850 ``'meas'``, rather than ``'ref'``.
852 Parameters
853 ----------
854 col : `str`
855 Name of flux column from which to compute magnitude.
856 Can be parseable by the `~lsst.pipe.tasks.functors.fluxName` function;
857 that is, you can pass ``'modelfit_CModel'`` instead of
858 ``'modelfit_CModel_instFlux'``, and it will understand.
859 """
860 _defaultDataset = 'meas'
862 def __init__(self, col, **kwargs):
863 self.col = fluxName(col)
864 # TO DO: DM-21955 Replace hard coded photometic calibration values.
865 self.fluxMag0 = 63095734448.0194
867 super().__init__(**kwargs)
869 @property
870 def columns(self):
871 return [self.col]
873 def _func(self, df):
874 with warnings.catch_warnings():
875 warnings.filterwarnings('ignore', r'invalid value encountered')
876 warnings.filterwarnings('ignore', r'divide by zero')
877 return -2.5*np.log10(df[self.col] / self.fluxMag0)
879 @property
880 def name(self):
881 return f'mag_{self.col}'
884class MagErr(Mag):
885 """Compute calibrated magnitude uncertainty.
887 Parameters
888 ----------
889 col : `str`
890 Name of the flux column.
891 """
893 def __init__(self, *args, **kwargs):
894 super().__init__(*args, **kwargs)
895 # TO DO: DM-21955 Replace hard coded photometic calibration values.
896 self.fluxMag0Err = 0.
898 @property
899 def columns(self):
900 return [self.col, self.col + 'Err']
902 def _func(self, df):
903 with warnings.catch_warnings():
904 warnings.filterwarnings('ignore', r'invalid value encountered')
905 warnings.filterwarnings('ignore', r'divide by zero')
906 fluxCol, fluxErrCol = self.columns
907 x = df[fluxErrCol] / df[fluxCol]
908 y = self.fluxMag0Err / self.fluxMag0
909 magErr = (2.5 / np.log(10.)) * np.sqrt(x*x + y*y)
910 return magErr
912 @property
913 def name(self):
914 return super().name + '_err'
917class MagDiff(Functor):
918 """Functor to calculate magnitude difference."""
919 _defaultDataset = 'meas'
921 def __init__(self, col1, col2, **kwargs):
922 self.col1 = fluxName(col1)
923 self.col2 = fluxName(col2)
924 super().__init__(**kwargs)
926 @property
927 def columns(self):
928 return [self.col1, self.col2]
930 def _func(self, df):
931 with warnings.catch_warnings():
932 warnings.filterwarnings('ignore', r'invalid value encountered')
933 warnings.filterwarnings('ignore', r'divide by zero')
934 return -2.5*np.log10(df[self.col1]/df[self.col2])
936 @property
937 def name(self):
938 return f'(mag_{self.col1} - mag_{self.col2})'
940 @property
941 def shortname(self):
942 return f'magDiff_{self.col1}_{self.col2}'
945class Color(Functor):
946 """Compute the color between two filters.
948 Computes color by initializing two different `Mag` functors based on the
949 ``col`` and filters provided, and then returning the difference.
951 This is enabled by the `_func` method expecting a DataFrame with a
952 multilevel column index, with both ``'band'`` and ``'column'``, instead of
953 just ``'column'``, which is the `Functor` default.
954 This is controlled by the `_dfLevels` attribute.
956 Also of note, the default dataset for `Color` is ``forced_src'``, whereas
957 for `Mag` it is ``'meas'``.
959 Parameters
960 ----------
961 col : str
962 Name of the flux column from which to compute; same as would be passed
963 to `~lsst.pipe.tasks.functors.Mag`.
965 filt2, filt1 : str
966 Filters from which to compute magnitude difference.
967 Color computed is ``Mag(filt2) - Mag(filt1)``.
968 """
969 _defaultDataset = 'forced_src'
970 _dfLevels = ('band', 'column')
971 _defaultNoDup = True
973 def __init__(self, col, filt2, filt1, **kwargs):
974 self.col = fluxName(col)
975 if filt2 == filt1: 975 ↛ 976line 975 didn't jump to line 976 because the condition on line 975 was never true
976 raise RuntimeError("Cannot compute Color for %s: %s - %s " % (col, filt2, filt1))
977 self.filt2 = filt2
978 self.filt1 = filt1
980 self.mag2 = Mag(col, filt=filt2, **kwargs)
981 self.mag1 = Mag(col, filt=filt1, **kwargs)
983 super().__init__(**kwargs)
985 @property
986 def filt(self):
987 return None
989 @filt.setter
990 def filt(self, filt):
991 pass
993 def _func(self, df):
994 mag2 = self.mag2._func(df[self.filt2])
995 mag1 = self.mag1._func(df[self.filt1])
996 return mag2 - mag1
998 @property
999 def columns(self):
1000 return [self.mag1.col, self.mag2.col]
1002 def multilevelColumns(self, parq, **kwargs):
1003 return [(self.dataset, self.filt1, self.col), (self.dataset, self.filt2, self.col)]
1005 @property
1006 def name(self):
1007 return f'{self.filt2} - {self.filt1} ({self.col})'
1009 @property
1010 def shortname(self):
1011 return f"{self.col}_{self.filt2.replace('-', '')}m{self.filt1.replace('-', '')}"
1014class DeconvolvedMoments(Functor):
1015 """This functor subtracts the trace of the PSF second moments from the
1016 trace of the second moments of the source.
1018 If the HsmShapeAlgorithm measurement is valid, then these will be used for
1019 the sources.
1020 Otherwise, the SdssShapeAlgorithm measurements will be used.
1021 """
1022 name = 'Deconvolved Moments'
1023 shortname = 'deconvolvedMoments'
1024 _columns = ("ext_shapeHSM_HsmSourceMoments_xx",
1025 "ext_shapeHSM_HsmSourceMoments_yy",
1026 "base_SdssShape_xx", "base_SdssShape_yy",
1027 "ext_shapeHSM_HsmPsfMoments_xx",
1028 "ext_shapeHSM_HsmPsfMoments_yy")
1030 def _func(self, df):
1031 """Calculate deconvolved moments."""
1032 if "ext_shapeHSM_HsmSourceMoments_xx" in df.columns: # _xx added by tdm 1032 ↛ 1035line 1032 didn't jump to line 1035 because the condition on line 1032 was always true
1033 hsm = df["ext_shapeHSM_HsmSourceMoments_xx"] + df["ext_shapeHSM_HsmSourceMoments_yy"]
1034 else:
1035 hsm = np.ones(len(df))*np.nan
1036 sdss = df["base_SdssShape_xx"] + df["base_SdssShape_yy"]
1037 if "ext_shapeHSM_HsmPsfMoments_xx" in df.columns: 1037 ↛ 1043line 1037 didn't jump to line 1043 because the condition on line 1037 was always true
1038 psf = df["ext_shapeHSM_HsmPsfMoments_xx"] + df["ext_shapeHSM_HsmPsfMoments_yy"]
1039 else:
1040 # LSST does not have shape.sdss.psf.
1041 # We could instead add base_PsfShape to the catalog using
1042 # exposure.getPsf().computeShape(s.getCentroid()).getIxx().
1043 raise RuntimeError('No psf shape parameter found in catalog')
1045 return hsm.where(np.isfinite(hsm), sdss) - psf
1048class SdssTraceSize(Functor):
1049 """Functor to calculate the SDSS trace radius size for sources.
1051 The SDSS trace radius size is a measure of size equal to the square root of
1052 half of the trace of the second moments tensor measured with the
1053 SdssShapeAlgorithm plugin.
1054 This has units of pixels.
1055 """
1056 name = "SDSS Trace Size"
1057 shortname = 'sdssTrace'
1058 _columns = ("base_SdssShape_xx", "base_SdssShape_yy")
1060 def _func(self, df):
1061 srcSize = np.sqrt(0.5*(df["base_SdssShape_xx"] + df["base_SdssShape_yy"]))
1062 return srcSize
1065class PsfSdssTraceSizeDiff(Functor):
1066 """Functor to calculate the SDSS trace radius size difference (%) between
1067 the object and the PSF model.
1069 See Also
1070 --------
1071 SdssTraceSize
1072 """
1073 name = "PSF - SDSS Trace Size"
1074 shortname = 'psf_sdssTrace'
1075 _columns = ("base_SdssShape_xx", "base_SdssShape_yy",
1076 "base_SdssShape_psf_xx", "base_SdssShape_psf_yy")
1078 def _func(self, df):
1079 srcSize = np.sqrt(0.5*(df["base_SdssShape_xx"] + df["base_SdssShape_yy"]))
1080 psfSize = np.sqrt(0.5*(df["base_SdssShape_psf_xx"] + df["base_SdssShape_psf_yy"]))
1081 sizeDiff = 100*(srcSize - psfSize)/(0.5*(srcSize + psfSize))
1082 return sizeDiff
1085class HsmTraceSize(Functor):
1086 """Functor to calculate the HSM trace radius size for sources.
1088 The HSM trace radius size is a measure of size equal to the square root of
1089 half of the trace of the second moments tensor measured with the
1090 HsmShapeAlgorithm plugin.
1091 This has units of pixels.
1092 """
1093 name = 'HSM Trace Size'
1094 shortname = 'hsmTrace'
1095 _columns = ("ext_shapeHSM_HsmSourceMoments_xx",
1096 "ext_shapeHSM_HsmSourceMoments_yy")
1098 def _func(self, df):
1099 srcSize = np.sqrt(0.5*(df["ext_shapeHSM_HsmSourceMoments_xx"]
1100 + df["ext_shapeHSM_HsmSourceMoments_yy"]))
1101 return srcSize
1104class PsfHsmTraceSizeDiff(Functor):
1105 """Functor to calculate the HSM trace radius size difference (%) between
1106 the object and the PSF model.
1108 See Also
1109 --------
1110 HsmTraceSize
1111 """
1112 name = 'PSF - HSM Trace Size'
1113 shortname = 'psf_HsmTrace'
1114 _columns = ("ext_shapeHSM_HsmSourceMoments_xx",
1115 "ext_shapeHSM_HsmSourceMoments_yy",
1116 "ext_shapeHSM_HsmPsfMoments_xx",
1117 "ext_shapeHSM_HsmPsfMoments_yy")
1119 def _func(self, df):
1120 srcSize = np.sqrt(0.5*(df["ext_shapeHSM_HsmSourceMoments_xx"]
1121 + df["ext_shapeHSM_HsmSourceMoments_yy"]))
1122 psfSize = np.sqrt(0.5*(df["ext_shapeHSM_HsmPsfMoments_xx"]
1123 + df["ext_shapeHSM_HsmPsfMoments_yy"]))
1124 sizeDiff = 100*(srcSize - psfSize)/(0.5*(srcSize + psfSize))
1125 return sizeDiff
1128class HsmFwhm(Functor):
1129 """Functor to calculate the PSF FWHM with second moments measured from the
1130 HsmShapeAlgorithm plugin.
1132 This is in units of arcseconds, and assumes the hsc_rings_v1 skymap pixel
1133 scale of 0.168 arcseconds/pixel.
1135 Notes
1136 -----
1137 This conversion assumes the PSF is Gaussian, which is not always the case.
1138 """
1139 name = 'HSM Psf FWHM'
1140 _columns = ('ext_shapeHSM_HsmPsfMoments_xx', 'ext_shapeHSM_HsmPsfMoments_yy')
1141 # TODO: DM-21403 pixel scale should be computed from the CD matrix or transform matrix
1142 pixelScale = 0.168
1143 SIGMA2FWHM = 2*np.sqrt(2*np.log(2))
1145 def _func(self, df):
1146 return (self.pixelScale*self.SIGMA2FWHM*np.sqrt(
1147 0.5*(df['ext_shapeHSM_HsmPsfMoments_xx']
1148 + df['ext_shapeHSM_HsmPsfMoments_yy']))).astype(np.float32)
1151class E1(Functor):
1152 r"""Calculate :math:`e_1` ellipticity component for sources, defined as:
1154 .. math::
1155 e_1 &= (I_{xx}-I_{yy})/(I_{xx}+I_{yy})
1157 See Also
1158 --------
1159 E2
1160 """
1161 name = "Distortion Ellipticity (e1)"
1162 shortname = "Distortion"
1164 def __init__(self, colXX, colXY, colYY, **kwargs):
1165 self.colXX = colXX
1166 self.colXY = colXY
1167 self.colYY = colYY
1168 self._columns = [self.colXX, self.colXY, self.colYY]
1169 super().__init__(**kwargs)
1171 @property
1172 def columns(self):
1173 return [self.colXX, self.colXY, self.colYY]
1175 def _func(self, df):
1176 return ((df[self.colXX] - df[self.colYY]) / (
1177 df[self.colXX] + df[self.colYY])).astype(np.float32)
1180class E2(Functor):
1181 r"""Calculate :math:`e_2` ellipticity component for sources, defined as:
1183 .. math::
1184 e_2 &= 2I_{xy}/(I_{xx}+I_{yy})
1186 See Also
1187 --------
1188 E1
1189 """
1190 name = "Ellipticity e2"
1192 def __init__(self, colXX, colXY, colYY, **kwargs):
1193 self.colXX = colXX
1194 self.colXY = colXY
1195 self.colYY = colYY
1196 super().__init__(**kwargs)
1198 @property
1199 def columns(self):
1200 return [self.colXX, self.colXY, self.colYY]
1202 def _func(self, df):
1203 return (2*df[self.colXY] / (df[self.colXX] + df[self.colYY])).astype(np.float32)
1206class RadiusFromQuadrupole(Functor):
1207 """Calculate the radius from the quadrupole moments.
1209 This returns the fourth root of the determinant of the second moments
1210 tensor, which has units of pixels.
1212 See Also
1213 --------
1214 SdssTraceSize
1215 HsmTraceSize
1216 """
1218 def __init__(self, colXX, colXY, colYY, **kwargs):
1219 self.colXX = colXX
1220 self.colXY = colXY
1221 self.colYY = colYY
1222 super().__init__(**kwargs)
1224 @property
1225 def columns(self):
1226 return [self.colXX, self.colXY, self.colYY]
1228 def _func(self, df):
1229 return ((df[self.colXX]*df[self.colYY] - df[self.colXY]**2)**0.25).astype(np.float32)
1232class LocalWcs(Functor):
1233 """Computations using the stored localWcs."""
1234 name = "LocalWcsOperations"
1236 def __init__(self,
1237 colCD_1_1,
1238 colCD_1_2,
1239 colCD_2_1,
1240 colCD_2_2,
1241 **kwargs):
1242 self.colCD_1_1 = colCD_1_1
1243 self.colCD_1_2 = colCD_1_2
1244 self.colCD_2_1 = colCD_2_1
1245 self.colCD_2_2 = colCD_2_2
1246 super().__init__(**kwargs)
1248 def computeDeltaRaDec(self, x, y, cd11, cd12, cd21, cd22):
1249 """Compute the dRA, dDec from dx, dy.
1251 Parameters
1252 ----------
1253 x : `~pandas.Series`
1254 X pixel coordinate.
1255 y : `~pandas.Series`
1256 Y pixel coordinate.
1257 cd11 : `~pandas.Series`
1258 [1, 1] element of the local Wcs affine transform.
1259 cd12 : `~pandas.Series`
1260 [1, 2] element of the local Wcs affine transform.
1261 cd21 : `~pandas.Series`
1262 [2, 1] element of the local Wcs affine transform.
1263 cd22 : `~pandas.Series`
1264 [2, 2] element of the local Wcs affine transform.
1266 Returns
1267 -------
1268 raDecTuple : tuple
1269 RA and Dec conversion of x and y given the local Wcs.
1270 Returned units are in radians.
1272 Notes
1273 -----
1274 If x and y are with respect to the CRVAL1, CRVAL2
1275 then this will return the RA, Dec for that WCS.
1276 """
1277 return (x * cd11 + y * cd12, x * cd21 + y * cd22)
1279 def computeSkySeparation(self, ra1, dec1, ra2, dec2):
1280 """Compute the local pixel scale conversion.
1282 Parameters
1283 ----------
1284 ra1 : `~pandas.Series`
1285 Ra of the first coordinate in radians.
1286 dec1 : `~pandas.Series`
1287 Dec of the first coordinate in radians.
1288 ra2 : `~pandas.Series`
1289 Ra of the second coordinate in radians.
1290 dec2 : `~pandas.Series`
1291 Dec of the second coordinate in radians.
1293 Returns
1294 -------
1295 dist : `~pandas.Series`
1296 Distance on the sphere in radians.
1297 """
1298 deltaDec = dec2 - dec1
1299 deltaRa = ra2 - ra1
1300 return 2 * np.arcsin(
1301 np.sqrt(
1302 np.sin(deltaDec / 2) ** 2
1303 + np.cos(dec2) * np.cos(dec1) * np.sin(deltaRa / 2) ** 2))
1305 def getSkySeparationFromPixel(self, x1, y1, x2, y2, cd11, cd12, cd21, cd22):
1306 """Compute the distance on the sphere from x2, y1 to x1, y1.
1308 Parameters
1309 ----------
1310 x1 : `~pandas.Series`
1311 X pixel coordinate.
1312 y1 : `~pandas.Series`
1313 Y pixel coordinate.
1314 x2 : `~pandas.Series`
1315 X pixel coordinate.
1316 y2 : `~pandas.Series`
1317 Y pixel coordinate.
1318 cd11 : `~pandas.Series`
1319 [1, 1] element of the local Wcs affine transform.
1320 cd12 : `~pandas.Series`
1321 [1, 2] element of the local Wcs affine transform.
1322 cd21 : `~pandas.Series`
1323 [2, 1] element of the local Wcs affine transform.
1324 cd22 : `~pandas.Series`
1325 [2, 2] element of the local Wcs affine transform.
1327 Returns
1328 -------
1329 Distance : `~pandas.Series`
1330 Arcseconds per pixel at the location of the local WC.
1331 """
1332 ra1, dec1 = self.computeDeltaRaDec(x1, y1, cd11, cd12, cd21, cd22)
1333 ra2, dec2 = self.computeDeltaRaDec(x2, y2, cd11, cd12, cd21, cd22)
1334 # Great circle distance for small separations.
1335 return self.computeSkySeparation(ra1, dec1, ra2, dec2)
1337 def computePositionAngle(self, ra1, dec1, ra2, dec2):
1338 """Compute position angle (E of N) from (ra1, dec1) to (ra2, dec2).
1340 Parameters
1341 ----------
1342 ra1 : iterable [`float`]
1343 RA of the first coordinate [radian].
1344 dec1 : iterable [`float`]
1345 Dec of the first coordinate [radian].
1346 ra2 : iterable [`float`]
1347 RA of the second coordinate [radian].
1348 dec2 : iterable [`float`]
1349 Dec of the second coordinate [radian].
1351 Returns
1352 -------
1353 Position Angle: `~pandas.Series`
1354 radians E of N
1356 Notes
1357 -----
1358 (ra1, dec1) -> (ra2, dec2) is interpreted as the shorter way around the sphere
1360 For a separation of 0.0001 rad, the position angle is good to 0.0009 rad
1361 all over the sphere.
1362 """
1363 # lsst.geom.SpherePoint has "bearingTo", which returns angle N of E
1364 # We instead want the astronomy convention of "Position Angle", which is angle E of N
1365 position_angle = np.zeros(len(ra1))
1366 for i, (r1, d1, r2, d2) in enumerate(zip(ra1, dec1, ra2, dec2)):
1367 point1 = geom.SpherePoint(r1, d1, geom.radians)
1368 point2 = geom.SpherePoint(r2, d2, geom.radians)
1369 bearing = point1.bearingTo(point2)
1370 pa_ref_angle = geom.Angle(np.pi/2, geom.radians) # in bearing system
1371 pa = pa_ref_angle - bearing
1372 # Wrap around to get Delta_RA from -pi to +pi
1373 pa = pa.wrapCtr()
1374 position_angle[i] = pa.asRadians()
1376 return pd.Series(position_angle)
1378 def getPositionAngleFromDetectorAngle(self, theta, cd11, cd12, cd21, cd22):
1379 """Compute position angle (E of N) from detector angle (+y of +x).
1381 Parameters
1382 ----------
1383 theta : `float`
1384 detector angle [radian]
1385 cd11 : `float`
1386 [1, 1] element of the local Wcs affine transform.
1387 cd12 : `float`
1388 [1, 2] element of the local Wcs affine transform.
1389 cd21 : `float`
1390 [2, 1] element of the local Wcs affine transform.
1391 cd22 : `float`
1392 [2, 2] element of the local Wcs affine transform.
1394 Returns
1395 -------
1396 Position Angle: `~pandas.Series`
1397 Degrees E of N.
1398 """
1399 # Create a unit vector in (x, y) along da
1400 dx = np.cos(theta)
1401 dy = np.sin(theta)
1402 ra1, dec1 = self.computeDeltaRaDec(0, 0, cd11, cd12, cd21, cd22)
1403 ra2, dec2 = self.computeDeltaRaDec(dx, dy, cd11, cd12, cd21, cd22)
1404 # Position angle of vector from (RA1, Dec1) to (RA2, Dec2)
1405 return np.rad2deg(self.computePositionAngle(ra1, dec1, ra2, dec2))
1408class ComputePixelScale(LocalWcs):
1409 """Compute the local pixel scale from the stored CDMatrix.
1410 """
1411 name = "PixelScale"
1413 @property
1414 def columns(self):
1415 return [self.colCD_1_1,
1416 self.colCD_1_2,
1417 self.colCD_2_1,
1418 self.colCD_2_2]
1420 def pixelScaleArcseconds(self, cd11, cd12, cd21, cd22):
1421 """Compute the local pixel to scale conversion in arcseconds.
1423 Parameters
1424 ----------
1425 cd11 : `~pandas.Series`
1426 [1, 1] element of the local Wcs affine transform in radians.
1427 cd11 : `~pandas.Series`
1428 [1, 1] element of the local Wcs affine transform in radians.
1429 cd12 : `~pandas.Series`
1430 [1, 2] element of the local Wcs affine transform in radians.
1431 cd21 : `~pandas.Series`
1432 [2, 1] element of the local Wcs affine transform in radians.
1433 cd22 : `~pandas.Series`
1434 [2, 2] element of the local Wcs affine transform in radians.
1436 Returns
1437 -------
1438 pixScale : `~pandas.Series`
1439 Arcseconds per pixel at the location of the local WC.
1440 """
1441 return 3600 * np.degrees(np.sqrt(np.fabs(cd11 * cd22 - cd12 * cd21)))
1443 def _func(self, df):
1444 return self.pixelScaleArcseconds(df[self.colCD_1_1],
1445 df[self.colCD_1_2],
1446 df[self.colCD_2_1],
1447 df[self.colCD_2_2])
1450class ConvertPixelToArcseconds(ComputePixelScale):
1451 """Convert a value in units of pixels to units of arcseconds."""
1453 def __init__(self,
1454 col,
1455 colCD_1_1,
1456 colCD_1_2,
1457 colCD_2_1,
1458 colCD_2_2,
1459 **kwargs):
1460 self.col = col
1461 super().__init__(colCD_1_1,
1462 colCD_1_2,
1463 colCD_2_1,
1464 colCD_2_2,
1465 **kwargs)
1467 @property
1468 def name(self):
1469 return f"{self.col}_asArcseconds"
1471 @property
1472 def columns(self):
1473 return [self.col,
1474 self.colCD_1_1,
1475 self.colCD_1_2,
1476 self.colCD_2_1,
1477 self.colCD_2_2]
1479 def _func(self, df):
1480 return df[self.col] * self.pixelScaleArcseconds(df[self.colCD_1_1],
1481 df[self.colCD_1_2],
1482 df[self.colCD_2_1],
1483 df[self.colCD_2_2])
1486class ConvertPixelSqToArcsecondsSq(ComputePixelScale):
1487 """Convert a value in units of pixels squared to units of arcseconds
1488 squared.
1489 """
1491 def __init__(self,
1492 col,
1493 colCD_1_1,
1494 colCD_1_2,
1495 colCD_2_1,
1496 colCD_2_2,
1497 **kwargs):
1498 self.col = col
1499 super().__init__(colCD_1_1,
1500 colCD_1_2,
1501 colCD_2_1,
1502 colCD_2_2,
1503 **kwargs)
1505 @property
1506 def name(self):
1507 return f"{self.col}_asArcsecondsSq"
1509 @property
1510 def columns(self):
1511 return [self.col,
1512 self.colCD_1_1,
1513 self.colCD_1_2,
1514 self.colCD_2_1,
1515 self.colCD_2_2]
1517 def _func(self, df):
1518 pixScale = self.pixelScaleArcseconds(df[self.colCD_1_1],
1519 df[self.colCD_1_2],
1520 df[self.colCD_2_1],
1521 df[self.colCD_2_2])
1522 return df[self.col] * pixScale * pixScale
1525class ConvertDetectorAngleToPositionAngle(LocalWcs):
1526 """Compute a position angle from a detector angle and the stored CDMatrix.
1528 Returns
1529 -------
1530 position angle : degrees
1531 """
1533 name = "PositionAngle"
1535 def __init__(
1536 self,
1537 theta_col,
1538 colCD_1_1,
1539 colCD_1_2,
1540 colCD_2_1,
1541 colCD_2_2,
1542 **kwargs
1543 ):
1544 self.theta_col = theta_col
1545 super().__init__(colCD_1_1, colCD_1_2, colCD_2_1, colCD_2_2, **kwargs)
1547 @property
1548 def columns(self):
1549 return [
1550 self.theta_col,
1551 self.colCD_1_1,
1552 self.colCD_1_2,
1553 self.colCD_2_1,
1554 self.colCD_2_2
1555 ]
1557 def _func(self, df):
1558 return self.getPositionAngleFromDetectorAngle(
1559 df[self.theta_col],
1560 df[self.colCD_1_1],
1561 df[self.colCD_1_2],
1562 df[self.colCD_2_1],
1563 df[self.colCD_2_2]
1564 )
1567class ConvertDetectorAngleErrToPositionAngleErr(Functor):
1568 """Convert a detector angle error to a position angle error.
1570 Returns
1571 -------
1572 position angle error : degrees
1573 """
1575 name = "PositionAngleErr"
1577 def __init__(self, theta_err_col, **kwargs):
1578 self.theta_err_col = theta_err_col
1579 super().__init__(**kwargs)
1581 @property
1582 def columns(self):
1583 return [self.theta_err_col]
1585 def _func(self, df):
1586 return np.rad2deg(df[self.theta_err_col])
1589class ReferenceBand(Functor):
1590 """Return the band used to seed multiband forced photometry.
1592 This functor is to be used on Object tables.
1593 It converts the boolean merge_measurements_{band} columns into a single
1594 string representing the first band for which merge_measurements_{band}
1595 is True.
1597 Assumes the default priority order of i, r, z, y, g, u.
1598 """
1599 name = 'Reference Band'
1600 shortname = 'refBand'
1602 band_order = ("i", "r", "z", "y", "g", "u")
1604 @property
1605 def columns(self):
1606 # Build the actual input column list, not hardcoded ugrizy
1607 bands = [band for band in self.band_order if band in self.bands]
1608 # In the unlikely scenario that users attempt to add non-ugrizy bands
1609 bands += [band for band in self.bands if band not in self.band_order]
1610 return [f"merge_measurement_{band}" for band in bands]
1612 def _func(self, df: pd.DataFrame) -> pd.Series:
1613 def getFilterAliasName(row):
1614 # Get column name with the max value (True > False).
1615 colName = row.idxmax()
1616 return colName.replace('merge_measurement_', '')
1618 # Skip columns that are unavailable, because this functor requests the
1619 # superset of bands that could be included in the object table.
1620 columns = [col for col in self.columns if col in df.columns]
1621 # Makes a Series of dtype object if df is empty.
1622 return df[columns].apply(getFilterAliasName, axis=1,
1623 result_type='reduce').astype('object')
1625 def __init__(self, bands: tuple[str] | list[str] | None = None, **kwargs):
1626 super().__init__(**kwargs)
1627 self.bands = self.band_order if bands is None else tuple(bands)
1630class Photometry(Functor):
1631 """Base class for Object table calibrated fluxes and magnitudes."""
1632 # AB to NanoJansky (3631 Jansky).
1633 AB_FLUX_SCALE = (0 * u.ABmag).to_value(u.nJy)
1634 LOG_AB_FLUX_SCALE = 12.56
1635 FIVE_OVER_2LOG10 = 1.085736204758129569
1636 # TO DO: DM-21955 Replace hard coded photometic calibration values.
1637 COADD_ZP = 27
1639 def __init__(self, colFlux, colFluxErr=None, **kwargs):
1640 self.vhypot = np.vectorize(self.hypot)
1641 self.col = colFlux
1642 self.colFluxErr = colFluxErr
1644 self.fluxMag0 = 1./np.power(10, -0.4*self.COADD_ZP)
1645 self.fluxMag0Err = 0.
1647 super().__init__(**kwargs)
1649 @property
1650 def columns(self):
1651 return [self.col]
1653 @property
1654 def name(self):
1655 return f'mag_{self.col}'
1657 @classmethod
1658 def hypot(cls, a, b):
1659 """Compute sqrt(a^2 + b^2) without under/overflow."""
1660 if np.abs(a) < np.abs(b):
1661 a, b = b, a
1662 if a == 0.:
1663 return 0.
1664 q = b/a
1665 return np.abs(a) * np.sqrt(1. + q*q)
1667 def dn2flux(self, dn, fluxMag0):
1668 """Convert instrumental flux to nanojanskys."""
1669 return (self.AB_FLUX_SCALE * dn / fluxMag0).astype(np.float32)
1671 def dn2mag(self, dn, fluxMag0):
1672 """Convert instrumental flux to AB magnitude."""
1673 with warnings.catch_warnings():
1674 warnings.filterwarnings('ignore', r'invalid value encountered')
1675 warnings.filterwarnings('ignore', r'divide by zero')
1676 return (-2.5 * np.log10(dn/fluxMag0)).astype(np.float32)
1678 def dn2fluxErr(self, dn, dnErr, fluxMag0, fluxMag0Err):
1679 """Convert instrumental flux error to nanojanskys."""
1680 retVal = self.vhypot(dn * fluxMag0Err, dnErr * fluxMag0)
1681 retVal *= self.AB_FLUX_SCALE / fluxMag0 / fluxMag0
1682 return retVal.astype(np.float32)
1684 def dn2MagErr(self, dn, dnErr, fluxMag0, fluxMag0Err):
1685 """Convert instrumental flux error to AB magnitude error."""
1686 retVal = self.dn2fluxErr(dn, dnErr, fluxMag0, fluxMag0Err) / self.dn2flux(dn, fluxMag0)
1687 return (self.FIVE_OVER_2LOG10 * retVal).astype(np.float32)
1690class NanoJansky(Photometry):
1691 """Convert instrumental flux to nanojanskys."""
1692 def _func(self, df):
1693 return self.dn2flux(df[self.col], self.fluxMag0)
1696class NanoJanskyErr(Photometry):
1697 """Convert instrumental flux error to nanojanskys."""
1698 @property
1699 def columns(self):
1700 return [self.col, self.colFluxErr]
1702 def _func(self, df):
1703 retArr = self.dn2fluxErr(df[self.col], df[self.colFluxErr], self.fluxMag0, self.fluxMag0Err)
1704 return pd.Series(retArr, index=df.index)
1707class LocalPhotometry(Functor):
1708 """Base class for calibrating the specified instrument flux column using
1709 the local photometric calibration.
1711 Parameters
1712 ----------
1713 instFluxCol : `str`
1714 Name of the instrument flux column.
1715 instFluxErrCol : `str`
1716 Name of the assocated error columns for ``instFluxCol``.
1717 photoCalibCol : `str`
1718 Name of local calibration column.
1719 photoCalibErrCol : `str`, optional
1720 Error associated with ``photoCalibCol``. Ignored and deprecated; will
1721 be removed after v29.
1723 See Also
1724 --------
1725 LocalNanojansky
1726 LocalNanojanskyErr
1727 """
1728 logNJanskyToAB = (1 * u.nJy).to_value(u.ABmag)
1730 def __init__(self,
1731 instFluxCol,
1732 instFluxErrCol,
1733 photoCalibCol,
1734 photoCalibErrCol=None,
1735 **kwargs):
1736 self.instFluxCol = instFluxCol
1737 self.instFluxErrCol = instFluxErrCol
1738 self.photoCalibCol = photoCalibCol
1739 # TODO[DM-49400]: remove this check and the argument it corresponds to.
1740 if photoCalibErrCol is not None: 1740 ↛ 1741line 1740 didn't jump to line 1741 because the condition on line 1740 was never true
1741 warnings.warn("The photoCalibErrCol argument is deprecated and will be removed after v29.",
1742 category=FutureWarning)
1743 super().__init__(**kwargs)
1745 def instFluxToNanojansky(self, instFlux, localCalib):
1746 """Convert instrument flux to nanojanskys.
1748 Parameters
1749 ----------
1750 instFlux : `~numpy.ndarray` or `~pandas.Series`
1751 Array of instrument flux measurements.
1752 localCalib : `~numpy.ndarray` or `~pandas.Series`
1753 Array of local photometric calibration estimates.
1755 Returns
1756 -------
1757 calibFlux : `~numpy.ndarray` or `~pandas.Series`
1758 Array of calibrated flux measurements.
1759 """
1760 return instFlux * localCalib
1762 def instFluxErrToNanojanskyErr(self, instFlux, instFluxErr, localCalib, localCalibErr=None):
1763 """Convert instrument flux to nanojanskys.
1765 Parameters
1766 ----------
1767 instFlux : `~numpy.ndarray` or `~pandas.Series`
1768 Array of instrument flux measurements. Ignored (accepted for
1769 backwards compatibility and consistency with magnitude-error
1770 calculation methods).
1771 instFluxErr : `~numpy.ndarray` or `~pandas.Series`
1772 Errors on associated ``instFlux`` values.
1773 localCalib : `~numpy.ndarray` or `~pandas.Series`
1774 Array of local photometric calibration estimates.
1775 localCalibErr : `~numpy.ndarray` or `~pandas.Series`, optional
1776 Errors on associated ``localCalib`` values. Ignored and deprecated;
1777 will be removed after v29.
1779 Returns
1780 -------
1781 calibFluxErr : `~numpy.ndarray` or `~pandas.Series`
1782 Errors on calibrated flux measurements.
1783 """
1784 # TODO[DM-49400]: remove this check and the argument it corresponds to.
1785 if localCalibErr is not None: 1785 ↛ 1786line 1785 didn't jump to line 1786 because the condition on line 1785 was never true
1786 warnings.warn("The localCalibErr argument is deprecated and will be removed after v29.",
1787 category=FutureWarning)
1788 return instFluxErr * localCalib
1790 def instFluxToMagnitude(self, instFlux, localCalib):
1791 """Convert instrument flux to nanojanskys.
1793 Parameters
1794 ----------
1795 instFlux : `~numpy.ndarray` or `~pandas.Series`
1796 Array of instrument flux measurements.
1797 localCalib : `~numpy.ndarray` or `~pandas.Series`
1798 Array of local photometric calibration estimates.
1800 Returns
1801 -------
1802 calibMag : `~numpy.ndarray` or `~pandas.Series`
1803 Array of calibrated AB magnitudes.
1804 """
1805 return -2.5 * np.log10(self.instFluxToNanojansky(instFlux, localCalib)) + self.logNJanskyToAB
1807 def instFluxErrToMagnitudeErr(self, instFlux, instFluxErr, localCalib, localCalibErr=None):
1808 """Convert instrument flux err to nanojanskys.
1810 Parameters
1811 ----------
1812 instFlux : `~numpy.ndarray` or `~pandas.Series`
1813 Array of instrument flux measurements.
1814 instFluxErr : `~numpy.ndarray` or `~pandas.Series`
1815 Errors on associated ``instFlux`` values.
1816 localCalib : `~numpy.ndarray` or `~pandas.Series`
1817 Array of local photometric calibration estimates.
1818 localCalibErr : `~numpy.ndarray` or `~pandas.Series`, optional
1819 Errors on associated ``localCalib`` values. Ignored and deprecated;
1820 will be removed after v29.
1822 Returns
1823 -------
1824 calibMagErr: `~numpy.ndarray` or `~pandas.Series`
1825 Error on calibrated AB magnitudes.
1826 """
1827 # TODO[DM-49400]: remove this check and the argument it corresponds to.
1828 if localCalibErr is not None: 1828 ↛ 1829line 1828 didn't jump to line 1829 because the condition on line 1828 was never true
1829 warnings.warn("The localCalibErr argument is deprecated and will be removed after v29.",
1830 category=FutureWarning)
1831 err = self.instFluxErrToNanojanskyErr(instFlux, instFluxErr, localCalib)
1832 return 2.5 / np.log(10) * err / self.instFluxToNanojansky(instFlux, instFluxErr)
1835class LocalNanojansky(LocalPhotometry):
1836 """Compute calibrated fluxes using the local calibration value.
1838 This returns units of nanojanskys.
1839 """
1841 @property
1842 def columns(self):
1843 return [self.instFluxCol, self.photoCalibCol]
1845 @property
1846 def name(self):
1847 return f'flux_{self.instFluxCol}'
1849 def _func(self, df):
1850 return self.instFluxToNanojansky(df[self.instFluxCol],
1851 df[self.photoCalibCol]).astype(np.float32)
1854class LocalNanojanskyErr(LocalPhotometry):
1855 """Compute calibrated flux errors using the local calibration value.
1857 This returns units of nanojanskys.
1858 """
1860 @property
1861 def columns(self):
1862 return [self.instFluxCol, self.instFluxErrCol, self.photoCalibCol]
1864 @property
1865 def name(self):
1866 return f'fluxErr_{self.instFluxCol}'
1868 def _func(self, df):
1869 return self.instFluxErrToNanojanskyErr(df[self.instFluxCol], df[self.instFluxErrCol],
1870 df[self.photoCalibCol]).astype(np.float32)
1873class LocalDipoleMeanFlux(LocalPhotometry):
1874 """Compute absolute mean of dipole fluxes.
1876 See Also
1877 --------
1878 LocalNanojansky
1879 LocalNanojanskyErr
1880 LocalDipoleMeanFluxErr
1881 LocalDipoleDiffFlux
1882 LocalDipoleDiffFluxErr
1883 """
1884 def __init__(self,
1885 instFluxPosCol,
1886 instFluxNegCol,
1887 instFluxPosErrCol,
1888 instFluxNegErrCol,
1889 photoCalibCol,
1890 # TODO[DM-49400]: remove this option; it's already deprecated (in super).
1891 photoCalibErrCol=None,
1892 **kwargs):
1893 self.instFluxNegCol = instFluxNegCol
1894 self.instFluxPosCol = instFluxPosCol
1895 self.instFluxNegErrCol = instFluxNegErrCol
1896 self.instFluxPosErrCol = instFluxPosErrCol
1897 self.photoCalibCol = photoCalibCol
1898 super().__init__(instFluxNegCol,
1899 instFluxNegErrCol,
1900 photoCalibCol,
1901 photoCalibErrCol,
1902 **kwargs)
1904 @property
1905 def columns(self):
1906 return [self.instFluxPosCol,
1907 self.instFluxNegCol,
1908 self.photoCalibCol]
1910 @property
1911 def name(self):
1912 return f'dipMeanFlux_{self.instFluxPosCol}_{self.instFluxNegCol}'
1914 def _func(self, df):
1915 return 0.5*(np.fabs(self.instFluxToNanojansky(df[self.instFluxNegCol], df[self.photoCalibCol]))
1916 + np.fabs(self.instFluxToNanojansky(df[self.instFluxPosCol], df[self.photoCalibCol])))
1919class LocalDipoleMeanFluxErr(LocalDipoleMeanFlux):
1920 """Compute the error on the absolute mean of dipole fluxes.
1922 See Also
1923 --------
1924 LocalNanojansky
1925 LocalNanojanskyErr
1926 LocalDipoleMeanFlux
1927 LocalDipoleDiffFlux
1928 LocalDipoleDiffFluxErr
1929 """
1931 @property
1932 def columns(self):
1933 return [self.instFluxPosCol,
1934 self.instFluxNegCol,
1935 self.instFluxPosErrCol,
1936 self.instFluxNegErrCol,
1937 self.photoCalibCol]
1939 @property
1940 def name(self):
1941 return f'dipMeanFluxErr_{self.instFluxPosCol}_{self.instFluxNegCol}'
1943 def _func(self, df):
1944 return 0.5*np.hypot(df[self.instFluxNegErrCol], df[self.instFluxPosErrCol]) * df[self.photoCalibCol]
1947class LocalDipoleDiffFlux(LocalDipoleMeanFlux):
1948 """Compute the absolute difference of dipole fluxes.
1950 Calculated value is (abs(pos) - abs(neg)).
1952 See Also
1953 --------
1954 LocalNanojansky
1955 LocalNanojanskyErr
1956 LocalDipoleMeanFlux
1957 LocalDipoleMeanFluxErr
1958 LocalDipoleDiffFluxErr
1959 """
1961 @property
1962 def columns(self):
1963 return [self.instFluxPosCol,
1964 self.instFluxNegCol,
1965 self.photoCalibCol]
1967 @property
1968 def name(self):
1969 return f'dipDiffFlux_{self.instFluxPosCol}_{self.instFluxNegCol}'
1971 def _func(self, df):
1972 return (np.fabs(self.instFluxToNanojansky(df[self.instFluxPosCol], df[self.photoCalibCol]))
1973 - np.fabs(self.instFluxToNanojansky(df[self.instFluxNegCol], df[self.photoCalibCol])))
1976class LocalDipoleDiffFluxErr(LocalDipoleMeanFlux):
1977 """Compute the error on the absolute difference of dipole fluxes.
1979 See Also
1980 --------
1981 LocalNanojansky
1982 LocalNanojanskyErr
1983 LocalDipoleMeanFlux
1984 LocalDipoleMeanFluxErr
1985 LocalDipoleDiffFlux
1986 """
1988 @property
1989 def columns(self):
1990 return [self.instFluxPosCol,
1991 self.instFluxNegCol,
1992 self.instFluxPosErrCol,
1993 self.instFluxNegErrCol,
1994 self.photoCalibCol]
1996 @property
1997 def name(self):
1998 return f'dipDiffFluxErr_{self.instFluxPosCol}_{self.instFluxNegCol}'
2000 def _func(self, df):
2001 return np.hypot(df[self.instFluxPosErrCol], df[self.instFluxNegErrCol]) * df[self.photoCalibCol]
2004class Ebv(Functor):
2005 """Compute E(B-V) from dustmaps.sfd."""
2006 _defaultDataset = 'ref'
2007 name = "E(B-V)"
2008 shortname = "ebv"
2010 def __init__(self, **kwargs):
2011 # Import is only needed for Ebv.
2012 # Suppress unnecessary .dustmapsrc log message on import.
2013 with open(os.devnull, "w") as devnull:
2014 with redirect_stdout(devnull):
2015 from dustmaps.sfd import SFDQuery
2016 self._columns = ['coord_ra', 'coord_dec']
2017 self.sfd = SFDQuery()
2018 super().__init__(**kwargs)
2020 def _func(self, df):
2021 coords = SkyCoord(df['coord_ra'].values * u.rad, df['coord_dec'].values * u.rad)
2022 ebv = self.sfd(coords)
2023 return pd.Series(ebv, index=df.index).astype('float32')
2026class MomentsBase(Functor):
2027 """Base class for functors that use shape moments and localWCS
2029 Attributes
2030 ----------
2031 is_covariance : bool
2032 Whether the shape columns are terms of a covariance matrix. If False,
2033 they will be assumed to be terms of a correlation matrix instead.
2034 """
2036 is_covariance: bool = True
2038 def __init__(self,
2039 shape_1_1,
2040 shape_2_2,
2041 shape_1_2,
2042 colCD_1_1,
2043 colCD_1_2,
2044 colCD_2_1,
2045 colCD_2_2,
2046 **kwargs):
2047 self.shape_1_1 = shape_1_1
2048 self.shape_2_2 = shape_2_2
2049 self.shape_1_2 = shape_1_2
2050 self.colCD_1_1 = colCD_1_1
2051 self.colCD_1_2 = colCD_1_2
2052 self.colCD_2_1 = colCD_2_1
2053 self.colCD_2_2 = colCD_2_2
2054 super().__init__(**kwargs)
2056 @property
2057 def columns(self):
2058 return [
2059 self.shape_1_1,
2060 self.shape_2_2,
2061 self.shape_1_2,
2062 ] + self.columns_ref
2064 @property
2065 def columns_ref(self):
2066 """Return columns that are needed from the ref table."""
2067 return [
2068 self.colCD_1_1,
2069 self.colCD_1_2,
2070 self.colCD_2_1,
2071 self.colCD_2_2]
2073 def compute_ellipse_terms(self, df, sky: bool = True):
2074 r"""Return terms commonly used for ellipse parameterization conversions.
2076 Parameters
2077 ----------
2078 df
2079 The data frame.
2080 sky
2081 Whether to compute the terms in sky coordinates.
2082 If False, XX, YY and XY moments are used instead of
2083 UU, VV and UV.
2085 Returns
2086 -------
2087 xx_p_yy
2088 The sum of the diagonal terms of the covariance.
2089 xx_m_yy
2090 The difference of the diagonal terms of the covariance.
2091 t2
2092 A term similar to the discriminant of the quadratic formula.
2093 """
2094 xx = self.sky_uu(df) if sky else self.get_xx(df)
2095 yy = self.sky_vv(df) if sky else self.get_yy(df)
2096 xx_m_yy = xx - yy
2097 t2 = xx_m_yy**2 + 4.0*(self.sky_uv(df) if sky else self.get_xy(df))**2
2098 # TODO: Check alternative form that may be more stable for computing
2099 # the minor axis size (see gauss2d/src/ellipse.cc)
2100 # t2 = xx**2 + yy**2 - 2*(xx*yy - 2*xy**2)
2101 return xx + yy, xx_m_yy, t2
2103 def get_xx(self, df):
2104 xx = df[self.shape_1_1]
2105 return xx if self.is_covariance else xx**2
2107 def get_yy(self, df):
2108 yy = df[self.shape_2_2]
2109 return yy if self.is_covariance else yy**2
2111 def get_xy(self, df):
2112 xy = df[self.shape_1_2]
2113 return xy if self.is_covariance else xy*df[self.shape_1_1]*df[self.shape_2_2]
2115 # Each of sky_uu, sky_vv, sky_uv evaluates one element of
2116 # CD_matrix * moments_matrix * CD_matrix.T
2117 def sky_uu(self, df):
2118 """Return the component of the moments tensor aligned with the RA axis, in radians."""
2119 i_xx = self.get_xx(df)
2120 i_yy = self.get_yy(df)
2121 i_xy = self.get_xy(df)
2122 CD_1_1 = df[self.colCD_1_1]
2123 CD_1_2 = df[self.colCD_1_2]
2124 return (CD_1_1*(i_xx*CD_1_1 + i_xy*CD_1_2)
2125 + CD_1_2*(i_xy*CD_1_1 + i_yy*CD_1_2))
2127 def sky_vv(self, df):
2128 """Return the component of the moments tensor aligned with the dec axis, in radians."""
2129 i_xx = self.get_xx(df)
2130 i_yy = self.get_yy(df)
2131 i_xy = self.get_xy(df)
2132 CD_2_1 = df[self.colCD_2_1]
2133 CD_2_2 = df[self.colCD_2_2]
2134 return (CD_2_1*(i_xx*CD_2_1 + i_xy*CD_2_2)
2135 + CD_2_2*(i_xy*CD_2_1 + i_yy*CD_2_2))
2137 def sky_uv(self, df):
2138 """Return the covariance of the moments tensor in ra, dec coordinates, in radians."""
2139 i_xx = self.get_xx(df)
2140 i_yy = self.get_yy(df)
2141 i_xy = self.get_xy(df)
2142 CD_1_1 = df[self.colCD_1_1]
2143 CD_1_2 = df[self.colCD_1_2]
2144 CD_2_1 = df[self.colCD_2_1]
2145 CD_2_2 = df[self.colCD_2_2]
2146 return ((CD_1_1 * i_xx + CD_1_2 * i_xy) * CD_2_1
2147 + (CD_1_1 * i_xy + CD_1_2 * i_yy) * CD_2_2)
2149 def get_g1(self, df):
2150 """
2151 Calculate shear-type ellipticity parameter G1.
2152 """
2153 # TODO: Replace this with functionality from afwGeom, DM-54015
2154 sky_uu = self.sky_uu(df)
2155 sky_vv = self.sky_vv(df)
2156 sky_uv = self.sky_uv(df)
2157 denom = sky_uu + sky_vv + 2 * np.sqrt(sky_uu*sky_vv - sky_uv**2)
2158 return ((sky_uu - sky_vv) / denom).astype(np.float32)
2160 def get_g2(self, df):
2161 """
2162 Calculate shear-type ellipticity parameter G2.
2164 This has the opposite sign as sky_uv in order to maintain consistency with the HSM moments
2165 sign convention.
2166 """
2167 # TODO: Replace this with functionality from afwGeom, DM-54015
2168 sky_uu = self.sky_uu(df)
2169 sky_vv = self.sky_vv(df)
2170 sky_uv = self.sky_uv(df)
2171 denom = sky_uu + sky_vv + 2 * np.sqrt(sky_uu*sky_vv - sky_uv**2)
2172 return (-2*sky_uv / denom).astype(np.float32)
2174 def get_trace(self, df):
2175 sky_uu = self.sky_uu(df)
2176 sky_vv = self.sky_vv(df)
2177 return np.sqrt(0.5*(sky_uu + sky_vv)).astype(np.float32)
2180class MomentsG1Sky(MomentsBase):
2181 """Rotate pixel moments Ixx,Iyy,Ixy into RA/dec frame and G1/G2 reduced
2182 shear parameterization"""
2183 _defaultDataset = 'meas'
2184 name = "moments_g1"
2185 shortname = "moments_g1"
2187 def _func(self, df):
2188 sky_g1 = self.get_g1(df)
2190 return pd.Series(sky_g1.astype(np.float32), index=df.index)
2193class MomentsG2Sky(MomentsBase):
2194 """Rotate pixel moments Ixx,Iyy,Ixy into RA/dec frame and G1/G2 reduced
2195 shear parameterization"""
2196 _defaultDataset = 'meas'
2197 name = "moments_g2"
2198 shortname = "moments_g2"
2200 def _func(self, df):
2201 sky_g2 = self.get_g2(df)
2203 return pd.Series(sky_g2.astype(np.float32), index=df.index)
2206class MomentsTraceSky(MomentsBase):
2207 """Trace radius size in arcseconds from pixel moments Ixx,Iyy,Ixy
2209 The trace radius size is a measure of size equal to the square root of
2210 half of the trace of the second moments tensor.
2211 """
2212 _defaultDataset = 'meas'
2213 name = "moments_trace"
2214 shortname = "moments_trace"
2216 def _func(self, df):
2217 sky_trace_radians = self.get_trace(df)
2219 return pd.Series((sky_trace_radians*(180/np.pi)*3600).astype(np.float32), index=df.index)
2222class MomentsIuuSky(MomentsBase):
2223 """Rotate pixel moments Ixx,Iyy,Ixy into ra,dec frame and arcseconds"""
2224 _defaultDataset = 'meas'
2225 name = "moments_uu"
2226 shortname = "moments_uu"
2228 def _func(self, df):
2229 sky_uu_radians = self.sky_uu(df)
2231 return pd.Series((sky_uu_radians*((180/np.pi)*3600)**2).astype(np.float32), index=df.index)
2234class CorrelationIuuSky(MomentsIuuSky):
2235 """MomentsIuuSky but from sigma_x, sigma_y, rho correlation terms."""
2236 is_covariance = False
2239class MomentsIvvSky(MomentsBase):
2240 """Rotate pixel moments Ixx,Iyy,Ixy into ra,dec frame and arcseconds"""
2241 _defaultDataset = 'meas'
2242 name = "moments_vv"
2243 shortname = "moments_vv"
2245 def _func(self, df):
2246 sky_vv_radians = self.sky_vv(df)
2248 return pd.Series((sky_vv_radians*((180/np.pi)*3600)**2).astype(np.float32), index=df.index)
2251class CorrelationIvvSky(MomentsIvvSky):
2252 """MomentsIvvSky but from sigma_x, sigma_y, rho correlation terms."""
2253 is_covariance = False
2256class MomentsIuvSky(MomentsBase):
2257 """Rotate pixel moments Ixx,Iyy,Ixy into ra,dec frame and arcseconds"""
2258 _defaultDataset = 'meas'
2259 name = "moments_uv"
2260 shortname = "moments_uv"
2262 def _func(self, df):
2263 sky_uv_radians = self.sky_uv(df)
2265 return pd.Series((sky_uv_radians*((180/np.pi)*3600)**2).astype(np.float32), index=df.index)
2268class CorrelationIuvSky(MomentsIuvSky):
2269 """MomentsIuvSky but from sigma_x, sigma_y, rho correlation terms."""
2270 is_covariance = False
2273class PositionAngleFromMoments(MomentsBase):
2274 """Compute position angle relative to ra,dec frame, in degrees, from Ixx,Iyy,Ixy pixel moments."""
2275 _defaultDataset = 'meas'
2276 name = "moments_theta"
2277 shortname = "moments_theta"
2279 def _func(self, df):
2280 sky_uu = self.sky_uu(df)
2281 sky_vv = self.sky_vv(df)
2282 sky_uv = self.sky_uv(df)
2283 theta = 0.5*np.arctan2(2*sky_uv, sky_uu - sky_vv)
2285 return pd.Series((np.degrees(np.array(theta))).astype(np.float32), index=df.index)
2288class PositionAngleFromCorrelation(PositionAngleFromMoments):
2289 """PositionAngleFromMoments but from sigma_x, sigma_y, rho correlation terms."""
2290 is_covariance = False
2293class SemimajorAxisFromMoments(MomentsBase):
2294 """Compute the semimajor axis length in arcseconds, from Ixx,Iyy,Ixy pixel moments."""
2295 _defaultDataset = 'meas'
2296 name = "moments_a"
2297 shortname = "moments_a"
2299 def _func(self, df):
2300 xx_p_yy, _, t2 = self.compute_ellipse_terms(df)
2301 # This copies what is done (unvectorized) in afw.geom.ellipse
2302 a_radians = np.sqrt(0.5 * (xx_p_yy + np.sqrt(t2)))
2304 return pd.Series((np.degrees(a_radians)*3600).astype(np.float32), index=df.index)
2307class SemimajorAxisFromCorrelation(SemimajorAxisFromMoments):
2308 """SemimajorAxisFromMoments but from sigma_x, sigma_y, rho correlation terms."""
2309 is_covariance = False
2312class SemiminorAxisFromMoments(MomentsBase):
2313 """Compute the semiminor axis length in arcseconds, from Ixx,Iyy,Ixy pixel moments."""
2314 _defaultDataset = 'meas'
2315 name = "moments_b"
2316 shortname = "moments_b"
2318 def _func(self, df):
2319 xx_p_yy, _, t2 = self.compute_ellipse_terms(df)
2320 # This copies what is done (unvectorized) in afw.geom.ellipse
2321 b_radians = np.sqrt(0.5 * (xx_p_yy - np.sqrt(t2)))
2323 return pd.Series((np.degrees(b_radians)*3600).astype(np.float32), index=df.index)
2326class SemiminorAxisFromCorrelation(SemiminorAxisFromMoments):
2327 """SemiminorAxisFromMoments but from sigma_x, sigma_y, rho correlation terms."""
2328 is_covariance = False