Skip to content

Commit

Permalink
Add array-based operations to basic photocalib conversions.
Browse files Browse the repository at this point in the history
  • Loading branch information
erykoff committed Mar 10, 2022
1 parent 761dc73 commit 998eb60
Show file tree
Hide file tree
Showing 3 changed files with 144 additions and 2 deletions.
1 change: 1 addition & 0 deletions python/lsst/afw/image/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,5 +30,6 @@
from ._exposureSummaryStats import *
from .basicUtils import *
from .testUtils import *
from ._photoCalibContinued import *

from ._exposureFitsReaderContinued import * # just here to support deprecation
103 changes: 103 additions & 0 deletions python/lsst/afw/image/_photoCalibContinued.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,103 @@
# This file is part of afw.
#
# Developed for the LSST Data Management System.
# This product includes software developed by the LSST Project
# (https://www.lsst.org).
# See the COPYRIGHT file at the top-level directory of this distribution
# for details of code ownership.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.

__all__ = ["PhotoCalib"]

import numpy as np
from astropy import units

from lsst.utils import continueClass

from ._imageLib import PhotoCalib


@continueClass
class PhotoCalib: # noqa: F811
def getLocalCalibrationArray(self, x, y):
"""Get the local calibration values for numpy arrays.
Parameters
----------
x : `np.ndarray` (N,)
Array of x values.
y : `np.ndarray` (N,)
Array of y values.
Returns
-------
localCalibration : `np.ndarray` (N,)
Array of local calibration values.
"""
if self.isConstant:
return np.full(len(x), self.getCalibrationMean())
else:
bf = self.computeScaledCalibration()
return self.getCalibrationMean()*bf.evaluate(x, y)

def instFluxToMagnitudeArray(self, instFluxes, x, y):
"""Convert instFlux to magnitudes for numpy arrays.
Parameters
----------
instFluxes : `np.ndarray` (N,)
Array of instFluxes to convert.
x : `np.ndarray` (N,)
Array of x values.
y : `np.ndarray` (N,)
Array of y values.
Returns
-------
magnitudes : `astropy.units.Magnitude` (N,)
Array of AB magnitudes.
"""
scale = self.getLocalCalibrationArray(x, y)
nanoJansky = (instFluxes*scale)*units.nJy

return nanoJansky.to(units.ABmag)

def magnitudeToInstFluxArray(self, magnitudes, x, y):
"""Convert magnitudes to instFlux for numpy arrays.
Parameters
----------
magnitudes : `np.ndarray` or `astropy.units.Magnitude` (N,)
Array of AB magnitudes.
x : `np.ndarray` (N,)
Array of x values.
y : `np.ndarray` (N,)
Array of y values.
Returns
-------
instFluxes : `np.ndarray` (N,)
Array of instFluxes.
"""
scale = self.getLocalCalibrationArray(x, y)

if not isinstance(magnitudes, units.Magnitude):
_magnitudes = magnitudes*units.ABmag
else:
_magnitudes = magnitudes

nanoJansky = _magnitudes.to(units.nJy).value

return nanoJansky/scale
42 changes: 40 additions & 2 deletions tests/test_photoCalib.py
Original file line number Diff line number Diff line change
Expand Up @@ -253,17 +253,55 @@ def _testPhotoCalibCenter(self, photoCalib, calibrationErr):
self.assertFloatsAlmostEqual(self.instFlux1,
photoCalib.magnitudeToInstFlux(mag, self.pointXShift),
rtol=1e-15)
mag = photoCalib.instFluxToMagnitude(self.instFlux2, self.pointXShift)
mag2 = photoCalib.instFluxToMagnitude(self.instFlux2, self.pointXShift)
self.assertFloatsAlmostEqual(self.instFlux2,
photoCalib.magnitudeToInstFlux(mag, self.pointXShift),
photoCalib.magnitudeToInstFlux(mag2, self.pointXShift),
rtol=1e-15)

# test round-tripping arrays (position specified)
instFlux1Array = np.full(10, self.instFlux1)
instFlux2Array = np.full(10, self.instFlux2)
pointXShiftXArray = np.full(10, self.pointXShift.getX())
pointXShiftYArray = np.full(10, self.pointXShift.getY())

magArray = photoCalib.instFluxToMagnitudeArray(
instFlux1Array,
pointXShiftXArray,
pointXShiftYArray
)
self.assertFloatsAlmostEqual(magArray.value, mag)
self.assertFloatsAlmostEqual(photoCalib.magnitudeToInstFluxArray(magArray,
pointXShiftXArray,
pointXShiftYArray
),
instFlux1Array,
rtol=5e-15)
mag2Array = photoCalib.instFluxToMagnitudeArray(
np.full(10, self.instFlux2),
np.full(10, self.pointXShift.getX()),
np.full(10, self.pointXShift.getY())
)
self.assertFloatsAlmostEqual(mag2Array.value, mag2)
self.assertFloatsAlmostEqual(photoCalib.magnitudeToInstFluxArray(mag2Array,
pointXShiftXArray,
pointXShiftYArray
),
instFlux2Array,
rtol=5e-15)

# test getLocalCalibration.
meas = photoCalib.instFluxToNanojansky(self.instFlux1, self.instFluxErr1, self.pointXShift)
localCalib = photoCalib.getLocalCalibration(self.pointXShift)
flux = localCalib * self.instFlux1
self.assertAlmostEqual(meas.value, flux)

# test getLocalCalibrationArray
localCalib2 = photoCalib.getLocalCalibrationArray(
pointXShiftXArray,
pointXShiftYArray
)
np.testing.assert_array_equal(localCalib2, localCalib)

def _testSourceCatalog(self, photoCalib, catalog, expectNanojansky, expectMag):
"""Test instFluxTo*(sourceCatalog, ...), and calibrateCatalog()."""

Expand Down

0 comments on commit 998eb60

Please sign in to comment.