LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+1b3060144d,g18429d2f64+f642bf4753,g199a45376c+0ba108daf9,g1fd858c14a+2dcf163641,g262e1987ae+7b8c96d2ca,g29ae962dfc+3bd6ecb08a,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53e1a9e7c5,g4595892280+fef73a337f,g47891489e3+2efcf17695,g4d44eb3520+642b70b07e,g53246c7159+8c5ae1fdc5,g67b6fd64d1+2efcf17695,g67fd3c3899+b70e05ef52,g74acd417e5+317eb4c7d4,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+2efcf17695,g8d7436a09f+3be3c13596,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+a4e7b16d9b,g97be763408+ad77d7208f,g9dd6db0277+b70e05ef52,ga681d05dcb+a3f46e7fff,gabf8522325+735880ea63,gac2eed3f23+2efcf17695,gb89ab40317+2efcf17695,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+b70e05ef52,gdab6d2f7ff+317eb4c7d4,gdc713202bf+b70e05ef52,gdfd2d52018+b10e285e0f,ge365c994fd+310e8507c4,ge410e46f29+2efcf17695,geaed405ab2+562b3308c0,gffca2db377+8c5ae1fdc5,w.2025.35
LSST Data Management Base Package
Loading...
Searching...
No Matches
PhotoCalib.cc
Go to the documentation of this file.
1/*
2 * LSST Data Management System
3 * Copyright 2017 LSST Corporation.
4 *
5 * This product includes software developed by the
6 * LSST Project (http://www.lsst.org/).
7 *
8 * This program is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the LSST License Statement and
19 * the GNU General Public License along with this program. If not,
20 * see <http://www.lsstcorp.org/LegalNotices/>.
21 */
22
23#include <cmath>
24#include <iostream>
25#include <iomanip>
26
27#include "lsst/geom/Point.h"
33#include "lsst/pex/exceptions.h"
34#include "ndarray.h"
37
38namespace lsst {
39namespace afw {
40
41template std::shared_ptr<image::PhotoCalib> table::io::PersistableFacade<image::PhotoCalib>::dynamicCast(
42 std::shared_ptr<table::io::Persistable> const &);
43
44namespace image {
45
46// ------------------- helpers -------------------
47
51 s << "value=" << measurement.value << ", error=" << measurement.error;
52 return os << s.str();
53}
54
55namespace {
56
57int const SERIALIZATION_VERSION = 1;
58
59double toNanojansky(double instFlux, double scale) { return instFlux * scale; }
60
61double toMagnitude(double instFlux, double scale) {
62 return cpputils::nanojanskyToABMagnitude(instFlux * scale);
63}
64
65double toInstFluxFromMagnitude(double magnitude, double scale) {
66 // Note: flux[nJy] / scale = instFlux[counts]
67 return cpputils::ABMagnitudeToNanojansky(magnitude) / scale;
68}
69
70double toInstFluxFromNanojansky(double nanojansky, double scale) { return nanojansky / scale; }
71
72double toNanojanskyErr(double instFluxErr, double scale) {
73 return instFluxErr * scale;
74}
75
76double toMagnitudeErr(double instFlux, double instFluxErr) {
77 return 2.5 / std::log(10.0) * (instFluxErr / instFlux);
78}
79
80} // anonymous namespace
81
82// ------------------- Conversions to nanojansky -------------------
83
84double PhotoCalib::instFluxToNanojansky(double instFlux, lsst::geom::Point<double, 2> const &point) const {
85 return toNanojansky(instFlux, evaluate(point));
86}
87
88double PhotoCalib::instFluxToNanojansky(double instFlux) const {
89 return toNanojansky(instFlux, _calibrationMean);
90}
91
92Measurement PhotoCalib::instFluxToNanojansky(double instFlux, double instFluxErr,
93 lsst::geom::Point<double, 2> const &point) const {
94 double calibration, error, nanojansky;
95 calibration = evaluate(point);
96 nanojansky = toNanojansky(instFlux, calibration);
97 error = toNanojanskyErr(instFluxErr, calibration);
98 return Measurement(nanojansky, error);
99}
100
101Measurement PhotoCalib::instFluxToNanojansky(double instFlux, double instFluxErr) const {
102 double nanojansky = toNanojansky(instFlux, _calibrationMean);
103 double error = toNanojanskyErr(instFluxErr, _calibrationMean);
104 return Measurement(nanojansky, error);
105}
106
108 std::string const &instFluxField) const {
109 auto position = sourceRecord.getCentroid();
110 auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFlux").key;
111 auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFluxErr").key;
112 return instFluxToNanojansky(sourceRecord.get(instFluxKey), sourceRecord.get(instFluxErrKey), position);
113}
114ndarray::Array<double, 2, 2> PhotoCalib::instFluxToNanojansky(afw::table::SourceCatalog const &sourceCatalog,
115 std::string const &instFluxField) const {
116 ndarray::Array<double, 2, 2> result =
117 ndarray::allocate(ndarray::makeVector(int(sourceCatalog.size()), 2));
118 instFluxToNanojanskyArray(sourceCatalog, instFluxField, result);
119 return result;
120}
121
123 std::string const &instFluxField, std::string const &outField) const {
124 auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
125 auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
126 auto nanojanskyKey = sourceCatalog.getSchema().find<double>(outField + "_flux").key;
127 auto nanojanskyErrKey = sourceCatalog.getSchema().find<double>(outField + "_fluxErr").key;
128 for (auto &record : sourceCatalog) {
129 auto result = instFluxToNanojansky(record.get(instFluxKey), record.get(instFluxErrKey),
130 record.getCentroid());
131 record.set(nanojanskyKey, result.value);
132 record.set(nanojanskyErrKey, result.error);
133 }
134}
135
136// ------------------- Conversions to Magnitudes -------------------
137
138double PhotoCalib::instFluxToMagnitude(double instFlux, lsst::geom::Point<double, 2> const &point) const {
139 return toMagnitude(instFlux, evaluate(point));
140}
141
142double PhotoCalib::instFluxToMagnitude(double instFlux) const {
143 return toMagnitude(instFlux, _calibrationMean);
144}
145
146Measurement PhotoCalib::instFluxToMagnitude(double instFlux, double instFluxErr,
147 lsst::geom::Point<double, 2> const &point) const {
148 double calibration, error, magnitude;
149 calibration = evaluate(point);
150 magnitude = toMagnitude(instFlux, calibration);
151 error = toMagnitudeErr(instFlux, instFluxErr);
152 return Measurement(magnitude, error);
153}
154
155Measurement PhotoCalib::instFluxToMagnitude(double instFlux, double instFluxErr) const {
156 double magnitude = toMagnitude(instFlux, _calibrationMean);
157 double error = toMagnitudeErr(instFlux, instFluxErr);
158 return Measurement(magnitude, error);
159}
160
162 std::string const &instFluxField) const {
163 auto position = sourceRecord.getCentroid();
164 auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFlux").key;
165 auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFluxErr").key;
166 return instFluxToMagnitude(sourceRecord.get(instFluxKey), sourceRecord.get(instFluxErrKey), position);
167}
168
169ndarray::Array<double, 2, 2> PhotoCalib::instFluxToMagnitude(afw::table::SourceCatalog const &sourceCatalog,
170 std::string const &instFluxField) const {
171 ndarray::Array<double, 2, 2> result =
172 ndarray::allocate(ndarray::makeVector(int(sourceCatalog.size()), 2));
173 instFluxToMagnitudeArray(sourceCatalog, instFluxField, result);
174 return result;
175}
176
178 std::string const &instFluxField, std::string const &outField) const {
179 auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
180 auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
181 auto magKey = sourceCatalog.getSchema().find<double>(outField + "_mag").key;
182 auto magErrKey = sourceCatalog.getSchema().find<double>(outField + "_magErr").key;
183 for (auto &record : sourceCatalog) {
184 auto result = instFluxToMagnitude(record.get(instFluxKey), record.get(instFluxErrKey),
185 record.getCentroid());
186 record.set(magKey, result.value);
187 record.set(magErrKey, result.error);
188 }
189}
190
191// ------------------- other utility methods -------------------
192
193double PhotoCalib::magnitudeToInstFlux(double magnitude) const {
194 return toInstFluxFromMagnitude(magnitude, _calibrationMean);
195}
196
197double PhotoCalib::magnitudeToInstFlux(double magnitude, lsst::geom::Point<double, 2> const &point) const {
198 return toInstFluxFromMagnitude(magnitude, evaluate(point));
199}
200
201double PhotoCalib::nanojanskyToInstFlux(double nanojansky) const {
202 return toInstFluxFromNanojansky(nanojansky, _calibrationMean);
203}
204
205double PhotoCalib::nanojanskyToInstFlux(double nanojansky, lsst::geom::Point<double, 2> const &point) const {
206 return toInstFluxFromNanojansky(nanojansky, evaluate(point));
207}
208
210 return *(_calibration) / _calibrationMean;
211}
212
216
217bool PhotoCalib::operator==(PhotoCalib const &rhs) const {
218 return (_calibrationMean == rhs._calibrationMean && _calibrationErr == rhs._calibrationErr &&
219 (*_calibration) == *(rhs._calibration));
220}
221
222double PhotoCalib::computeCalibrationMean(std::shared_ptr<afw::math::BoundedField> calibration) const {
223 return calibration->mean();
224}
225
227 return std::make_unique<PhotoCalib>(*this);
228}
229
231 std::stringstream buffer;
232 if (_isConstant)
233 buffer << "spatially constant with ";
234 else
235 buffer << *_calibration << " with ";
236 buffer << "mean: " << _calibrationMean << " error: " << _calibrationErr;
237 return buffer.str();
238}
239
240bool PhotoCalib::equals(typehandling::Storable const &other) const noexcept {
241 return singleClassEquals(*this, other);
242}
243
245 return os << photoCalib.toString();
246}
247
249 // Deep copy construct, as we're multiplying in-place.
250 auto result = MaskedImage<float>(maskedImage, true);
251 if (_isConstant) {
252 result *= _calibrationMean;
253 } else {
254 _calibration->multiplyImage(result, true); // only in the overlap region
255 }
256 return result;
257}
258
259
261 bool includeScaleUncertainty) const {
262 return calibrateImage(maskedImage);
263}
264
266 // Deep copy construct, as we're multiplying in-place.
267 auto result = MaskedImage<float>(maskedImage, true);
268 if (_isConstant) {
269 result /= _calibrationMean;
270 } else {
271 _calibration->divideImage(result, true); // only in the overlap region
272 }
273 return result;
274}
275
277 bool includeScaleUncertainty) const {
278 return uncalibrateImage(maskedImage);
279}
280
282 std::vector<std::string> const &instFluxFields) const {
283 auto const &inSchema = catalog.getSchema();
284 afw::table::SchemaMapper mapper(inSchema, true); // true: share the alias map
285 mapper.addMinimalSchema(inSchema);
286
287 using FieldD = afw::table::Field<double>;
288
289 struct Keys {
290 table::Key<double> instFlux;
291 table::Key<double> instFluxErr;
293 table::Key<double> fluxErr;
295 table::Key<double> magErr;
296 };
297
299 keys.reserve(instFluxFields.size());
300 for (auto const &field : instFluxFields) {
301 Keys newKey;
302 newKey.instFlux = inSchema[inSchema.join(field, "instFlux")];
303 newKey.flux =
304 mapper.addOutputField(FieldD(inSchema.join(field, "flux"), "calibrated flux", "nJy"), true);
305 newKey.mag = mapper.addOutputField(
306 FieldD(inSchema.join(field, "mag"), "calibrated magnitude", "mag(AB)"), true);
307 try {
308 newKey.instFluxErr = inSchema.find<double>(inSchema.join(field, "instFluxErr")).key;
309 newKey.fluxErr = mapper.addOutputField(
310 FieldD(inSchema.join(field, "fluxErr"), "calibrated flux uncertainty", "nJy"), true);
311 newKey.magErr = mapper.addOutputField(
312 FieldD(inSchema.join(field, "magErr"), "calibrated magnitude uncertainty", "mag(AB)"),
313 true);
315 ; // Keys struct defaults to invalid keys; that marks the error as missing.
316 }
317 keys.emplace_back(newKey);
318 }
319
320 // Create the new catalog
322 output.insert(mapper, output.begin(), catalog.begin(), catalog.end());
323
324 auto calibration = evaluateCatalog(output);
325
326 // fill in the catalog values
327 int iRec = 0;
328 for (auto &rec : output) {
329 for (auto &key : keys) {
330 double instFlux = rec.get(key.instFlux);
331 double nanojansky = toNanojansky(instFlux, calibration[iRec]);
332 rec.set(key.flux, nanojansky);
333 rec.set(key.mag, toMagnitude(instFlux, calibration[iRec]));
334 if (key.instFluxErr.isValid()) {
335 double instFluxErr = rec.get(key.instFluxErr);
336 rec.set(key.fluxErr, toNanojanskyErr(instFluxErr, calibration[iRec]));
337 rec.set(key.magErr, toMagnitudeErr(instFlux, instFluxErr));
338 }
339 }
340 ++iRec;
341 }
342
343 return output;
344}
345
347 std::vector<std::string> instFluxFields;
348 static std::string const SUFFIX = "_instFlux";
349 for (auto const &name : catalog.getSchema().getNames()) {
350 // Pick every field ending in "_instFlux", grabbing everything before that prefix.
351 if (name.size() > SUFFIX.size() + 1 &&
352 name.compare(name.size() - SUFFIX.size(), SUFFIX.size(), SUFFIX) == 0) {
353 instFluxFields.emplace_back(name.substr(0, name.size() - 9));
354 }
355 }
356 return calibrateCatalog(catalog, instFluxFields);
357}
358
359// ------------------- persistence -------------------
360
361namespace {
362
363class PhotoCalibSchema {
364public:
365 table::Schema schema;
366 table::Key<double> calibrationMean;
367 table::Key<double> calibrationErr;
368 table::Key<table::Flag> isConstant;
369 table::Key<int> field;
371
372 // No copying
373 PhotoCalibSchema(PhotoCalibSchema const &) = delete;
374 PhotoCalibSchema &operator=(PhotoCalibSchema const &) = delete;
375 // No moving
376 PhotoCalibSchema(PhotoCalibSchema &&) = delete;
377 PhotoCalibSchema &operator=(PhotoCalibSchema &&) = delete;
378
379 static PhotoCalibSchema const &get() {
380 static PhotoCalibSchema const instance;
381 return instance;
382 }
383
384private:
385 PhotoCalibSchema()
386 : schema(),
387 calibrationMean(schema.addField<double>(
388 "calibrationMean", "mean calibration on this PhotoCalib's domain", "count")),
389 calibrationErr(
390 schema.addField<double>("calibrationErr", "1-sigma error on calibrationMean", "count")),
391 isConstant(schema.addField<table::Flag>("isConstant", "Is this spatially-constant?")),
392 field(schema.addField<int>("field", "archive ID of the BoundedField object")),
393 version(schema.addField<int>("version", "version of this PhotoCalib")) {}
394};
395
396class PhotoCalibFactory : public table::io::PersistableFactory {
397public:
398 std::shared_ptr<table::io::Persistable> read(InputArchive const &archive,
399 CatalogVector const &catalogs) const override {
400 table::BaseRecord const &record = catalogs.front().front();
401 PhotoCalibSchema const &keys = PhotoCalibSchema::get();
402 int version = getVersion(record);
403 if (version < 1) {
404 throw(pex::exceptions::RuntimeError("Unsupported version (version 0 was defined in maggies): " +
405 std::to_string(version)));
406 }
407 return std::make_shared<PhotoCalib>(record.get(keys.calibrationMean), record.get(keys.calibrationErr),
408 archive.get<afw::math::BoundedField>(record.get(keys.field)),
409 record.get(keys.isConstant));
410 }
411
412 PhotoCalibFactory(std::string const &name) : afw::table::io::PersistableFactory(name) {}
413
414protected:
415 int getVersion(table::BaseRecord const &record) const {
416 int version = -1;
417 try {
418 std::string versionName = "version";
419 auto versionKey = record.getSchema().find<int>(versionName);
420 version = record.get(versionKey.key);
421 } catch (const pex::exceptions::NotFoundError &) {
422 // un-versioned files are version 0
423 version = 0;
424 }
425 return version;
426 }
427};
428
429std::string getPhotoCalibPersistenceName() { return "PhotoCalib"; }
430
431PhotoCalibFactory registration(getPhotoCalibPersistenceName());
432
433} // namespace
434
435/*
436 * Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
437 */
438
439namespace {
440int const CALIB_TABLE_CURRENT_VERSION = 2; // final version of Calib in ExposureTable
441std::string const EXPTIME_FIELD_NAME = "exptime"; // name of exposure time field (no longer used)
442
443class CalibKeys {
444public:
447 table::Key<double> expTime;
448 table::Key<double> fluxMag0;
449 table::Key<double> fluxMag0Err;
450
451 // No copying
452 CalibKeys(const CalibKeys &) = delete;
453 CalibKeys &operator=(const CalibKeys &) = delete;
454
455 // No moving
456 CalibKeys(CalibKeys &&) = delete;
457 CalibKeys &operator=(CalibKeys &&) = delete;
458
459 CalibKeys(int tableVersion = CALIB_TABLE_CURRENT_VERSION)
460 : schema(), midTime(), expTime(), fluxMag0(), fluxMag0Err() {
461 if (tableVersion == 1) {
462 // obsolete fields
463 midTime = schema.addField<std::int64_t>(
464 "midtime", "middle of the time of the exposure relative to Unix epoch", "ns");
465 expTime = schema.addField<double>(EXPTIME_FIELD_NAME, "exposure time", "s");
466 }
467 fluxMag0 = schema.addField<double>("fluxmag0", "flux of a zero-magnitude object", "count");
468 fluxMag0Err = schema.addField<double>("fluxmag0.err", "1-sigma error on fluxmag0", "count");
469 }
470};
471
472class CalibFactory : public table::io::PersistableFactory {
473public:
474 std::shared_ptr<table::io::Persistable> read(InputArchive const &archive,
475 CatalogVector const &catalogs) const override {
476 // table version is not persisted, so we don't have a clean way to determine the version;
477 // the hack is version = 1 if exptime found, else current
478 int tableVersion = 1;
479 try {
480 catalogs.front().getSchema().find<double>(EXPTIME_FIELD_NAME);
481 } catch (pex::exceptions::NotFoundError const &) {
482 tableVersion = CALIB_TABLE_CURRENT_VERSION;
483 }
484
485 CalibKeys const keys{tableVersion};
486 LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
487 LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
488 LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
489 table::BaseRecord const &record = catalogs.front().front();
490
491 double calibration = cpputils::referenceFlux / record.get(keys.fluxMag0);
492 double calibrationErr = cpputils::referenceFlux * record.get(keys.fluxMag0Err) /
493 std::pow(record.get(keys.fluxMag0), 2);
494 return std::make_shared<PhotoCalib>(calibration, calibrationErr);
495 }
496
497 explicit CalibFactory(std::string const &name) : table::io::PersistableFactory(name) {}
498};
499
500std::string getCalibPersistenceName() { return "Calib"; }
501
502CalibFactory calibRegistration(getCalibPersistenceName());
503
504} // namespace
505
506std::string PhotoCalib::getPersistenceName() const { return getPhotoCalibPersistenceName(); }
507
509 PhotoCalibSchema const &keys = PhotoCalibSchema::get();
510 table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
511 auto record = catalog.addNew();
512 record->set(keys.calibrationMean, _calibrationMean);
513 record->set(keys.calibrationErr, _calibrationErr);
514 record->set(keys.isConstant, _isConstant);
515 record->set(keys.field, handle.put(_calibration));
516 record->set(keys.version, SERIALIZATION_VERSION);
517 handle.saveCatalog(catalog);
518}
519
520// ------------------- private/protected helpers -------------------
521
522double PhotoCalib::evaluate(lsst::geom::Point<double, 2> const &point) const {
523 if (_isConstant)
524 return _calibrationMean;
525 else
526 return _calibration->evaluate(point);
527}
528
529ndarray::Array<double, 1> PhotoCalib::evaluateArray(ndarray::Array<double, 1> const &xx,
530 ndarray::Array<double, 1> const &yy) const {
531 if (_isConstant) {
532 ndarray::Array<double, 1> result = ndarray::allocate(ndarray::makeVector(xx.size()));
533 result.deep() = _calibrationMean;
534 return result;
535 } else {
536 return _calibration->evaluate(xx, yy);
537 }
538}
539
540ndarray::Array<double, 1> PhotoCalib::evaluateCatalog(afw::table::SourceCatalog const &sourceCatalog) const {
541 ndarray::Array<double, 1> xx = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
542 ndarray::Array<double, 1> yy = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
543 size_t i = 0;
544 for (auto const &rec : sourceCatalog) {
545 auto point = rec.getCentroid();
546 xx[i] = point.getX();
547 yy[i] = point.getY();
548 ++i;
549 }
550 return evaluateArray(xx, yy);
551}
552
553void PhotoCalib::instFluxToNanojanskyArray(afw::table::SourceCatalog const &sourceCatalog,
554 std::string const &instFluxField,
555 ndarray::Array<double, 2, 2> result) const {
556 auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
557 auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
558
559 auto calibration = evaluateCatalog(sourceCatalog);
560 int i = 0;
561 auto iter = result.begin();
562 for (auto const &rec : sourceCatalog) {
563 double instFlux = rec.get(instFluxKey);
564 double instFluxErr = rec.get(instFluxErrKey);
565 double nanojansky = toNanojansky(instFlux, calibration[i]);
566 (*iter)[0] = nanojansky;
567 (*iter)[1] = toNanojanskyErr(instFluxErr, calibration[i]);
568 ++iter;
569 ++i;
570 }
571}
572
573void PhotoCalib::instFluxToMagnitudeArray(afw::table::SourceCatalog const &sourceCatalog,
574 std::string const &instFluxField,
575 ndarray::Array<double, 2, 2> result) const {
576 auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
577 auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
578
579 auto calibration = evaluateCatalog(sourceCatalog);
580 auto iter = result.begin();
581 int i = 0;
582 for (auto const &rec : sourceCatalog) {
583 double instFlux = rec.get(instFluxKey);
584 double instFluxErr = rec.get(instFluxErrKey);
585 (*iter)[0] = toMagnitude(instFlux, calibration[i]);
586 (*iter)[1] = toMagnitudeErr(instFlux, instFluxErr);
587 ++iter;
588 ++i;
589 }
590}
591
593 auto key = "FLUXMAG0";
594 if (metadata.exists(key)) {
595 double instFluxMag0 = metadata.getAsDouble(key);
596 if (strip) metadata.remove(key);
597
598 double instFluxMag0Err = 0.0;
599 key = "FLUXMAG0ERR";
600 if (metadata.exists(key)) {
601 instFluxMag0Err = metadata.getAsDouble(key);
602 if (strip) metadata.remove(key);
603 }
604 return makePhotoCalibFromCalibZeroPoint(instFluxMag0, instFluxMag0Err);
605 } else {
606 return nullptr;
607 }
608}
609
610std::shared_ptr<PhotoCalib> makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err) {
611 double calibration = cpputils::referenceFlux / instFluxMag0;
612 double calibrationErr = cpputils::referenceFlux * instFluxMag0Err / std::pow(instFluxMag0, 2);
613 return std::make_shared<PhotoCalib>(calibration, calibrationErr);
614}
615
616} // namespace image
617} // namespace afw
618} // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
Utilities for converting between flux and magnitude in C++.
Implementation of the Photometric Calibration class.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition Persistable.h:48
A class to manipulate images, masks, and variance as a single object.
Definition MaskedImage.h:74
The photometric calibration of an exposure.
Definition PhotoCalib.h:114
afw::table::SourceCatalog calibrateCatalog(afw::table::SourceCatalog const &catalog, std::vector< std::string > const &instFluxFields) const
Return a flux calibrated catalog, with new _flux, _fluxErr, _mag, and _magErr fields.
bool operator==(PhotoCalib const &rhs) const
Two PhotoCalibs are equal if their component bounded fields and calibrationErr are equal.
double instFluxToNanojansky(double instFlux, lsst::geom::Point< double, 2 > const &point) const
Convert instFlux in ADU to nJy at a point in the BoundedField.
Definition PhotoCalib.cc:84
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
PhotoCalib(PhotoCalib const &)=default
double instFluxToMagnitude(double instFlux, lsst::geom::Point< double, 2 > const &point) const
Convert instFlux in ADU to AB magnitude.
std::string toString() const override
Create a string representation of this object.
bool equals(typehandling::Storable const &other) const noexcept override
Compare this object to another Storable.
MaskedImage< float > uncalibrateImage(MaskedImage< float > const &maskedImage) const
Return a un-calibrated image, with pixel values in ADU (or whatever the original input to this photoC...
std::shared_ptr< typehandling::Storable > cloneStorable() const override
Create a new PhotoCalib that is a copy of this one.
double nanojanskyToInstFlux(double nanojansky, lsst::geom::Point< double, 2 > const &point) const
Convert nanojansky to instFlux (ADU).
std::shared_ptr< afw::math::BoundedField > computeScalingTo(std::shared_ptr< PhotoCalib > other) const
Calculates the scaling between this PhotoCalib and another PhotoCalib.
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
std::shared_ptr< afw::math::BoundedField > computeScaledCalibration() const
Calculates the spatially-variable calibration, normalized by the mean in the valid domain.
MaskedImage< float > calibrateImage(MaskedImage< float > const &maskedImage) const
Return a flux calibrated image, with pixel values in nJy.
double magnitudeToInstFlux(double magnitude, lsst::geom::Point< double, 2 > const &point) const
Convert AB magnitude to instFlux (ADU).
Base class for all records.
Definition BaseRecord.h:31
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition BaseRecord.h:151
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
Definition BaseRecord.h:80
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition Catalog.h:490
A class used as a handle to a particular field in a table.
Definition Key.h:53
Defines the fields and offsets for a table.
Definition Schema.h:51
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Definition Schema.cc:467
A mapping between the keys of two Schemas, used to copy data between them.
Schema const getOutputSchema() const
Return the output schema (copy-on-write).
Key< T > addOutputField(Field< T > const &newField, bool doReplace=false)
Add a new field to the output Schema that is not connected to the input Schema.
void addMinimalSchema(Schema const &minimal, bool doMap=true)
Add the given minimal schema to the output schema.
Record class that contains measurements made on a single exposure.
Definition Source.h:78
CentroidSlotDefinition::MeasValue getCentroid() const
Get the value of the Centroid slot measurement.
Definition Source.h:569
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
int put(Persistable const *obj, bool permissive=false)
Save an object to the archive and return a unique ID that can be used to retrieve it from an InputArc...
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
A base class for factory classes used to reconstruct objects from records.
io::OutputArchiveHandle OutputArchiveHandle
Interface supporting iteration over heterogenous containers.
Definition Storable.h:58
static bool singleClassEquals(T const &lhs, Storable const &rhs)
Test if a Storable is of a particular class and equal to another object.
Definition Storable.h:151
Class for storing generic metadata.
Definition PropertySet.h:67
virtual void remove(std::string const &name)
Remove all values for a property name (possibly hierarchical).
bool exists(std::string const &name) const
Determine if a name (possibly hierarchical) exists.
double getAsDouble(std::string const &name) const
Get the last value for any arithmetic property name (possibly hierarchical).
A coordinate class intended to represent absolute positions.
Definition Point.h:169
Reports errors in the logical structure of the program.
Definition Runtime.h:46
Reports attempts to access elements using an invalid key.
Definition Runtime.h:151
T emplace_back(T... args)
T get(T... args)
T log(T... args)
T make_shared(T... args)
scale(algorithm, min, max=None, frame=None)
Definition ds9.py:108
std::shared_ptr< PhotoCalib > makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err)
Construct a PhotoCalib from the deprecated Calib-style instFluxMag0/instFluxMag0Err values.
std::shared_ptr< PhotoCalib > makePhotoCalibFromMetadata(daf::base::PropertySet &metadata, bool strip=false)
Construct a PhotoCalib from FITS FLUXMAG0/FLUXMAG0ERR keywords.
std::ostream & operator<<(std::ostream &os, Measurement const &measurement)
Definition PhotoCalib.cc:48
CatalogT< BaseRecord > BaseCatalog
Definition fwd.h:72
SortedCatalogT< SourceRecord > SourceCatalog
Definition fwd.h:85
const double referenceFlux
The Oke & Gunn (1983) AB magnitude reference flux, in nJy (often approximated as 3631....
Definition Magnitude.h:46
double ABMagnitudeToNanojansky(double magnitude)
Convert an AB magnitude to a flux in nanojansky.
Definition Magnitude.cc:32
double nanojanskyToABMagnitude(double flux)
Convert a flux in nanojansky to AB magnitude.
Definition Magnitude.cc:30
T pow(T... args)
T setprecision(T... args)
T size(T... args)
T str(T... args)
A value and its error.
Definition PhotoCalib.h:51
A description of a field in a table.
Definition Field.h:24
T to_string(T... args)
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override