LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
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 
26 #include "lsst/geom/Point.h"
32 #include "lsst/pex/exceptions.h"
33 #include "ndarray.h"
35 #include "lsst/utils/Magnitude.h"
36 
37 namespace lsst {
38 namespace afw {
39 
42 
43 namespace image {
44 
45 // ------------------- helpers -------------------
46 
50  s << "value=" << measurement.value << ", error=" << measurement.error;
51  return os << s.str();
52 }
53 
54 namespace {
55 
56 int const SERIALIZATION_VERSION = 1;
57 
58 double toNanojansky(double instFlux, double scale) { return instFlux * scale; }
59 
60 double toMagnitude(double instFlux, double scale) { return utils::nanojanskyToABMagnitude(instFlux * scale); }
61 
62 double toInstFluxFromMagnitude(double magnitude, double scale) {
63  // Note: flux[nJy] / scale = instFlux[counts]
64  return utils::ABMagnitudeToNanojansky(magnitude) / scale;
65 }
66 
67 double toNanojanskyErr(double instFlux, double instFluxErr, double scale, double scaleErr,
68  double nanojansky) {
69  return std::abs(nanojansky) * hypot(instFluxErr / instFlux, scaleErr / scale);
70 }
71 
85 void toNanojanskyVariance(ndarray::Array<float const, 2, 1> const &instFlux,
86  ndarray::Array<float const, 2, 1> const &instFluxVar, float scaleErr,
87  ndarray::Array<float const, 2, 1> const &flux, ndarray::Array<float, 2, 1> out) {
88  auto eigenFlux = ndarray::asEigen<Eigen::ArrayXpr>(flux);
89  auto eigenInstFluxVar = ndarray::asEigen<Eigen::ArrayXpr>(instFluxVar);
90  auto eigenInstFlux = ndarray::asEigen<Eigen::ArrayXpr>(instFlux);
91  auto eigenOut = ndarray::asEigen<Eigen::ArrayXpr>(out);
92  eigenOut = eigenFlux.square() *
93  (eigenInstFluxVar / eigenInstFlux.square() + (scaleErr / eigenFlux * eigenInstFlux).square());
94 }
95 
96 double toMagnitudeErr(double instFlux, double instFluxErr, double scale, double scaleErr) {
97  return 2.5 / std::log(10.0) * hypot(instFluxErr / instFlux, scaleErr / scale);
98 }
99 
100 } // anonymous namespace
101 
102 // ------------------- Conversions to nanojansky -------------------
103 
104 double PhotoCalib::instFluxToNanojansky(double instFlux, lsst::geom::Point<double, 2> const &point) const {
105  return toNanojansky(instFlux, evaluate(point));
106 }
107 
108 double PhotoCalib::instFluxToNanojansky(double instFlux) const {
109  return toNanojansky(instFlux, _calibrationMean);
110 }
111 
112 Measurement PhotoCalib::instFluxToNanojansky(double instFlux, double instFluxErr,
113  lsst::geom::Point<double, 2> const &point) const {
114  double calibration, error, nanojansky;
115  calibration = evaluate(point);
116  nanojansky = toNanojansky(instFlux, calibration);
117  error = toNanojanskyErr(instFlux, instFluxErr, calibration, _calibrationErr, nanojansky);
118  return Measurement(nanojansky, error);
119 }
120 
121 Measurement PhotoCalib::instFluxToNanojansky(double instFlux, double instFluxErr) const {
122  double nanojansky = toNanojansky(instFlux, _calibrationMean);
123  double error = toNanojanskyErr(instFlux, instFluxErr, _calibrationMean, _calibrationErr, nanojansky);
124  return Measurement(nanojansky, error);
125 }
126 
127 Measurement PhotoCalib::instFluxToNanojansky(afw::table::SourceRecord const &sourceRecord,
128  std::string const &instFluxField) const {
129  auto position = sourceRecord.getCentroid();
130  auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFlux").key;
131  auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFluxErr").key;
132  return instFluxToNanojansky(sourceRecord.get(instFluxKey), sourceRecord.get(instFluxErrKey), position);
133 }
134 ndarray::Array<double, 2, 2> PhotoCalib::instFluxToNanojansky(afw::table::SourceCatalog const &sourceCatalog,
135  std::string const &instFluxField) const {
136  ndarray::Array<double, 2, 2> result =
137  ndarray::allocate(ndarray::makeVector(int(sourceCatalog.size()), 2));
138  instFluxToNanojanskyArray(sourceCatalog, instFluxField, result);
139  return result;
140 }
141 
142 void PhotoCalib::instFluxToNanojansky(afw::table::SourceCatalog &sourceCatalog,
143  std::string const &instFluxField, std::string const &outField) const {
144  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
145  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
146  auto nanojanskyKey = sourceCatalog.getSchema().find<double>(outField + "_flux").key;
147  auto nanojanskyErrKey = sourceCatalog.getSchema().find<double>(outField + "_fluxErr").key;
148  for (auto &record : sourceCatalog) {
149  auto result = instFluxToNanojansky(record.get(instFluxKey), record.get(instFluxErrKey),
150  record.getCentroid());
151  record.set(nanojanskyKey, result.value);
152  record.set(nanojanskyErrKey, result.error);
153  }
154 }
155 
156 // ------------------- Conversions to Magnitudes -------------------
157 
158 double PhotoCalib::instFluxToMagnitude(double instFlux, lsst::geom::Point<double, 2> const &point) const {
159  return toMagnitude(instFlux, evaluate(point));
160 }
161 
162 double PhotoCalib::instFluxToMagnitude(double instFlux) const {
163  return toMagnitude(instFlux, _calibrationMean);
164 }
165 
166 Measurement PhotoCalib::instFluxToMagnitude(double instFlux, double instFluxErr,
167  lsst::geom::Point<double, 2> const &point) const {
168  double calibration, error, magnitude;
169  calibration = evaluate(point);
170  magnitude = toMagnitude(instFlux, calibration);
171  error = toMagnitudeErr(instFlux, instFluxErr, calibration, _calibrationErr);
172  return Measurement(magnitude, error);
173 }
174 
175 Measurement PhotoCalib::instFluxToMagnitude(double instFlux, double instFluxErr) const {
176  double magnitude = toMagnitude(instFlux, _calibrationMean);
177  double error = toMagnitudeErr(instFlux, instFluxErr, _calibrationMean, _calibrationErr);
178  return Measurement(magnitude, error);
179 }
180 
181 Measurement PhotoCalib::instFluxToMagnitude(afw::table::SourceRecord const &sourceRecord,
182  std::string const &instFluxField) const {
183  auto position = sourceRecord.getCentroid();
184  auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFlux").key;
185  auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFluxErr").key;
186  return instFluxToMagnitude(sourceRecord.get(instFluxKey), sourceRecord.get(instFluxErrKey), position);
187 }
188 
189 ndarray::Array<double, 2, 2> PhotoCalib::instFluxToMagnitude(afw::table::SourceCatalog const &sourceCatalog,
190  std::string const &instFluxField) const {
191  ndarray::Array<double, 2, 2> result =
192  ndarray::allocate(ndarray::makeVector(int(sourceCatalog.size()), 2));
193  instFluxToMagnitudeArray(sourceCatalog, instFluxField, result);
194  return result;
195 }
196 
197 void PhotoCalib::instFluxToMagnitude(afw::table::SourceCatalog &sourceCatalog,
198  std::string const &instFluxField, std::string const &outField) const {
199  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
200  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
201  auto magKey = sourceCatalog.getSchema().find<double>(outField + "_mag").key;
202  auto magErrKey = sourceCatalog.getSchema().find<double>(outField + "_magErr").key;
203  for (auto &record : sourceCatalog) {
204  auto result = instFluxToMagnitude(record.get(instFluxKey), record.get(instFluxErrKey),
205  record.getCentroid());
206  record.set(magKey, result.value);
207  record.set(magErrKey, result.error);
208  }
209 }
210 
211 // ------------------- other utility methods -------------------
212 
213 double PhotoCalib::magnitudeToInstFlux(double magnitude) const {
214  return toInstFluxFromMagnitude(magnitude, _calibrationMean);
215 }
216 
217 double PhotoCalib::magnitudeToInstFlux(double magnitude, lsst::geom::Point<double, 2> const &point) const {
218  return toInstFluxFromMagnitude(magnitude, evaluate(point));
219 }
220 
221 std::shared_ptr<math::BoundedField> PhotoCalib::computeScaledCalibration() const {
222  return *(_calibration) / _calibrationMean;
223 }
224 
225 std::shared_ptr<math::BoundedField> PhotoCalib::computeScalingTo(std::shared_ptr<PhotoCalib> other) const {
226  throw LSST_EXCEPT(pex::exceptions::LogicError, "Not Implemented: See DM-10154.");
227 }
228 
229 bool PhotoCalib::operator==(PhotoCalib const &rhs) const {
230  return (_calibrationMean == rhs._calibrationMean && _calibrationErr == rhs._calibrationErr &&
231  (*_calibration) == *(rhs._calibration));
232 }
233 
234 double PhotoCalib::computeCalibrationMean(std::shared_ptr<afw::math::BoundedField> calibration) const {
235  return calibration->mean();
236 }
237 
238 std::shared_ptr<typehandling::Storable> PhotoCalib::cloneStorable() const {
239  return std::make_unique<PhotoCalib>(*this);
240 }
241 
242 std::string PhotoCalib::toString() const {
243  std::stringstream buffer;
244  if (_isConstant)
245  buffer << "spatially constant with ";
246  else
247  buffer << *_calibration << " with ";
248  buffer << "mean: " << _calibrationMean << " error: " << _calibrationErr;
249  return buffer.str();
250 }
251 
252 bool PhotoCalib::equals(typehandling::Storable const &other) const noexcept {
253  return singleClassEquals(*this, other);
254 }
255 
257  return os << photoCalib.toString();
258 }
259 
260 MaskedImage<float> PhotoCalib::calibrateImage(MaskedImage<float> const &maskedImage,
261  bool includeScaleUncertainty) const {
262  // Deep copy construct, as we're mutiplying in-place.
263  auto result = MaskedImage<float>(maskedImage, true);
264 
265  if (_isConstant) {
266  *(result.getImage()) *= _calibrationMean;
267  } else {
268  _calibration->multiplyImage(*(result.getImage()), true); // only in the overlap region
269  }
270  if (includeScaleUncertainty) {
271  toNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(),
272  _calibrationErr, result.getImage()->getArray(),
273  result.getVariance()->getArray());
274  } else {
275  toNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(), 0,
276  result.getImage()->getArray(), result.getVariance()->getArray());
277  }
278 
279  return result;
280 }
281 
282 afw::table::SourceCatalog PhotoCalib::calibrateCatalog(afw::table::SourceCatalog const &catalog,
283  std::vector<std::string> const &instFluxFields) const {
284  auto const &inSchema = catalog.getSchema();
285  afw::table::SchemaMapper mapper(inSchema, true); // true: share the alias map
286  mapper.addMinimalSchema(inSchema);
287 
288  using FieldD = afw::table::Field<double>;
289 
290  struct Keys {
291  table::Key<double> instFlux;
292  table::Key<double> instFluxErr;
293  table::Key<double> flux;
294  table::Key<double> fluxErr;
295  table::Key<double> mag;
296  table::Key<double> magErr;
297  };
298 
300  keys.reserve(instFluxFields.size());
301  for (auto const &field : instFluxFields) {
302  Keys newKey;
303  newKey.instFlux = inSchema[inSchema.join(field, "instFlux")];
304  newKey.flux =
305  mapper.addOutputField(FieldD(inSchema.join(field, "flux"), "calibrated flux", "nJy"), true);
306  newKey.mag = mapper.addOutputField(
307  FieldD(inSchema.join(field, "mag"), "calibrated magnitude", "mag(AB)"), true);
308  try {
309  newKey.instFluxErr = inSchema.find<double>(inSchema.join(field, "instFluxErr")).key;
310  newKey.fluxErr = mapper.addOutputField(
311  FieldD(inSchema.join(field, "fluxErr"), "calibrated flux uncertainty", "nJy"), true);
312  newKey.magErr = mapper.addOutputField(
313  FieldD(inSchema.join(field, "magErr"), "calibrated magnitude uncertainty", "mag(AB)"),
314  true);
315  } catch (pex::exceptions::NotFoundError &) {
316  ; // Keys struct defaults to invalid keys; that marks the error as missing.
317  }
318  keys.emplace_back(newKey);
319  }
320 
321  // Create the new catalog
322  afw::table::SourceCatalog output(mapper.getOutputSchema());
323  output.insert(mapper, output.begin(), catalog.begin(), catalog.end());
324 
325  auto calibration = evaluateCatalog(output);
326 
327  // fill in the catalog values
328  int iRec = 0;
329  for (auto &rec : output) {
330  for (auto &key : keys) {
331  double instFlux = rec.get(key.instFlux);
332  double nanojansky = toNanojansky(instFlux, calibration[iRec]);
333  rec.set(key.flux, nanojansky);
334  rec.set(key.mag, toMagnitude(instFlux, calibration[iRec]));
335  if (key.instFluxErr.isValid()) {
336  double instFluxErr = rec.get(key.instFluxErr);
337  rec.set(key.fluxErr, toNanojanskyErr(instFlux, instFluxErr, calibration[iRec],
338  _calibrationErr, nanojansky));
339  rec.set(key.magErr,
340  toMagnitudeErr(instFlux, instFluxErr, calibration[iRec], _calibrationErr));
341  }
342  }
343  ++iRec;
344  }
345 
346  return output;
347 }
348 
349 afw::table::SourceCatalog PhotoCalib::calibrateCatalog(afw::table::SourceCatalog const &catalog) const {
350  std::vector<std::string> instFluxFields;
351  static std::string const SUFFIX = "_instFlux";
352  for (auto const &name : catalog.getSchema().getNames()) {
353  // Pick every field ending in "_instFlux", grabbing everything before that prefix.
354  if (name.size() > SUFFIX.size() + 1 &&
355  name.compare(name.size() - SUFFIX.size(), SUFFIX.size(), SUFFIX) == 0) {
356  instFluxFields.emplace_back(name.substr(0, name.size() - 9));
357  }
358  }
359  return calibrateCatalog(catalog, instFluxFields);
360 }
361 
362 // ------------------- persistence -------------------
363 
364 namespace {
365 
366 class PhotoCalibSchema {
367 public:
368  table::Schema schema;
369  table::Key<double> calibrationMean;
370  table::Key<double> calibrationErr;
371  table::Key<table::Flag> isConstant;
372  table::Key<int> field;
373  table::Key<int> version;
374 
375  // No copying
376  PhotoCalibSchema(PhotoCalibSchema const &) = delete;
377  PhotoCalibSchema &operator=(PhotoCalibSchema const &) = delete;
378  // No moving
379  PhotoCalibSchema(PhotoCalibSchema &&) = delete;
380  PhotoCalibSchema &operator=(PhotoCalibSchema &&) = delete;
381 
382  static PhotoCalibSchema const &get() {
383  static PhotoCalibSchema const instance;
384  return instance;
385  }
386 
387 private:
388  PhotoCalibSchema()
389  : schema(),
390  calibrationMean(schema.addField<double>(
391  "calibrationMean", "mean calibration on this PhotoCalib's domain", "count")),
393  schema.addField<double>("calibrationErr", "1-sigma error on calibrationMean", "count")),
394  isConstant(schema.addField<table::Flag>("isConstant", "Is this spatially-constant?")),
395  field(schema.addField<int>("field", "archive ID of the BoundedField object")),
396  version(schema.addField<int>("version", "version of this PhotoCalib")) {}
397 };
398 
399 class PhotoCalibFactory : public table::io::PersistableFactory {
400 public:
402  read(InputArchive const &archive, CatalogVector const &catalogs) const override {
403  table::BaseRecord const &record = catalogs.front().front();
404  PhotoCalibSchema const &keys = PhotoCalibSchema::get();
405  int version = getVersion(record);
406  if (version < 1) {
407  throw(pex::exceptions::RuntimeError("Unsupported version (version 0 was defined in maggies): " +
409  }
410  return std::make_shared<PhotoCalib>(record.get(keys.calibrationMean), record.get(keys.calibrationErr),
411  archive.get<afw::math::BoundedField>(record.get(keys.field)),
412  record.get(keys.isConstant));
413  }
414 
415  PhotoCalibFactory(std::string const &name) : afw::table::io::PersistableFactory(name) {}
416 
417 protected:
418  int getVersion(table::BaseRecord const &record) const {
419  int version = -1;
420  try {
421  std::string versionName = "version";
422  auto versionKey = record.getSchema().find<int>(versionName);
423  version = record.get(versionKey.key);
424  } catch (const pex::exceptions::NotFoundError &) {
425  // un-versioned files are version 0
426  version = 0;
427  }
428  return version;
429  }
430 };
431 
432 std::string getPhotoCalibPersistenceName() { return "PhotoCalib"; }
433 
434 PhotoCalibFactory registration(getPhotoCalibPersistenceName());
435 
436 } // namespace
437 
442 namespace {
443 int const CALIB_TABLE_CURRENT_VERSION = 2; // final version of Calib in ExposureTable
444 std::string const EXPTIME_FIELD_NAME = "exptime"; // name of exposure time field (no longer used)
445 
446 class CalibKeys {
447 public:
448  table::Schema schema;
449  table::Key<std::int64_t> midTime;
450  table::Key<double> expTime;
451  table::Key<double> fluxMag0;
452  table::Key<double> fluxMag0Err;
453 
454  // No copying
455  CalibKeys(const CalibKeys &) = delete;
456  CalibKeys &operator=(const CalibKeys &) = delete;
457 
458  // No moving
459  CalibKeys(CalibKeys &&) = delete;
460  CalibKeys &operator=(CalibKeys &&) = delete;
461 
462  CalibKeys(int tableVersion = CALIB_TABLE_CURRENT_VERSION)
463  : schema(), midTime(), expTime(), fluxMag0(), fluxMag0Err() {
464  if (tableVersion == 1) {
465  // obsolete fields
467  "midtime", "middle of the time of the exposure relative to Unix epoch", "ns");
468  expTime = schema.addField<double>(EXPTIME_FIELD_NAME, "exposure time", "s");
469  }
470  fluxMag0 = schema.addField<double>("fluxmag0", "flux of a zero-magnitude object", "count");
471  fluxMag0Err = schema.addField<double>("fluxmag0.err", "1-sigma error on fluxmag0", "count");
472  }
473 };
474 
475 class CalibFactory : public table::io::PersistableFactory {
476 public:
477  std::shared_ptr<table::io::Persistable> read(InputArchive const &archive,
478  CatalogVector const &catalogs) const override {
479  // table version is not persisted, so we don't have a clean way to determine the version;
480  // the hack is version = 1 if exptime found, else current
481  int tableVersion = 1;
482  try {
483  catalogs.front().getSchema().find<double>(EXPTIME_FIELD_NAME);
484  } catch (pex::exceptions::NotFoundError const&) {
485  tableVersion = CALIB_TABLE_CURRENT_VERSION;
486  }
487 
488  CalibKeys const keys{tableVersion};
489  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
490  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
491  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
492  table::BaseRecord const &record = catalogs.front().front();
493 
494  double calibration = utils::referenceFlux / record.get(keys.fluxMag0);
495  double calibrationErr =
496  utils::referenceFlux * record.get(keys.fluxMag0Err) / std::pow(record.get(keys.fluxMag0), 2);
497  return std::make_shared<PhotoCalib>(calibration, calibrationErr);
498  }
499 
500  explicit CalibFactory(std::string const &name) : table::io::PersistableFactory(name) {}
501 };
502 
503 std::string getCalibPersistenceName() { return "Calib"; }
504 
505 CalibFactory calibRegistration(getCalibPersistenceName());
506 
507 } // namespace
508 
509 std::string PhotoCalib::getPersistenceName() const { return getPhotoCalibPersistenceName(); }
510 
512  PhotoCalibSchema const &keys = PhotoCalibSchema::get();
513  table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
514  auto record = catalog.addNew();
515  record->set(keys.calibrationMean, _calibrationMean);
516  record->set(keys.calibrationErr, _calibrationErr);
517  record->set(keys.isConstant, _isConstant);
518  record->set(keys.field, handle.put(_calibration));
519  record->set(keys.version, SERIALIZATION_VERSION);
520  handle.saveCatalog(catalog);
521 }
522 
523 // ------------------- private/protected helpers -------------------
524 
525 double PhotoCalib::evaluate(lsst::geom::Point<double, 2> const &point) const {
526  if (_isConstant)
527  return _calibrationMean;
528  else
529  return _calibration->evaluate(point);
530 }
531 
532 ndarray::Array<double, 1> PhotoCalib::evaluateArray(ndarray::Array<double, 1> const &xx,
533  ndarray::Array<double, 1> const &yy) const {
534  if (_isConstant) {
535  ndarray::Array<double, 1> result = ndarray::allocate(ndarray::makeVector(xx.size()));
536  result.deep() = _calibrationMean;
537  return result;
538  } else {
539  return _calibration->evaluate(xx, yy);
540  }
541 }
542 
543 ndarray::Array<double, 1> PhotoCalib::evaluateCatalog(afw::table::SourceCatalog const &sourceCatalog) const {
544  ndarray::Array<double, 1> xx = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
545  ndarray::Array<double, 1> yy = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
546  size_t i = 0;
547  for (auto const &rec : sourceCatalog) {
548  auto point = rec.getCentroid();
549  xx[i] = point.getX();
550  yy[i] = point.getY();
551  ++i;
552  }
553  return evaluateArray(xx, yy);
554 }
555 
556 void PhotoCalib::instFluxToNanojanskyArray(afw::table::SourceCatalog const &sourceCatalog,
557  std::string const &instFluxField,
558  ndarray::Array<double, 2, 2> result) const {
559  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
560  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
561 
562  auto calibration = evaluateCatalog(sourceCatalog);
563  int i = 0;
564  auto iter = result.begin();
565  for (auto const &rec : sourceCatalog) {
566  double instFlux = rec.get(instFluxKey);
567  double instFluxErr = rec.get(instFluxErrKey);
568  double nanojansky = toNanojansky(instFlux, calibration[i]);
569  (*iter)[0] = nanojansky;
570  (*iter)[1] = toNanojanskyErr(instFlux, instFluxErr, calibration[i], _calibrationErr, nanojansky);
571  ++iter;
572  ++i;
573  }
574 }
575 
576 void PhotoCalib::instFluxToMagnitudeArray(afw::table::SourceCatalog const &sourceCatalog,
577  std::string const &instFluxField,
578  ndarray::Array<double, 2, 2> result) const {
579  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
580  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
581 
582  auto calibration = evaluateCatalog(sourceCatalog);
583  auto iter = result.begin();
584  int i = 0;
585  for (auto const &rec : sourceCatalog) {
586  double instFlux = rec.get(instFluxKey);
587  double instFluxErr = rec.get(instFluxErrKey);
588  (*iter)[0] = toMagnitude(instFlux, calibration[i]);
589  (*iter)[1] = toMagnitudeErr(instFlux, instFluxErr, calibration[i], _calibrationErr);
590  ++iter;
591  ++i;
592  }
593 }
594 
596  auto key = "FLUXMAG0";
597  if (metadata.exists(key)) {
598  double instFluxMag0 = metadata.getAsDouble(key);
599  if (strip) metadata.remove(key);
600 
601  double instFluxMag0Err = 0.0;
602  key = "FLUXMAG0ERR";
603  if (metadata.exists(key)) {
604  instFluxMag0Err = metadata.getAsDouble(key);
605  if (strip) metadata.remove(key);
606  }
607  return makePhotoCalibFromCalibZeroPoint(instFluxMag0, instFluxMag0Err);
608  } else {
609  return nullptr;
610  }
611 }
612 
613 std::shared_ptr<PhotoCalib> makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err) {
614  double calibration = utils::referenceFlux / instFluxMag0;
615  double calibrationErr = utils::referenceFlux * instFluxMag0Err / std::pow(instFluxMag0, 2);
616  return std::make_shared<PhotoCalib>(calibration, calibrationErr);
617 }
618 
619 } // namespace image
620 } // namespace afw
621 } // namespace lsst
py::object result
Definition: _schema.cc:429
table::Key< std::string > name
Definition: Amplifier.cc:116
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
table::Key< double > calibrationMean
Definition: PhotoCalib.cc:369
table::Schema schema
Definition: PhotoCalib.cc:368
table::Key< std::int64_t > midTime
Definition: PhotoCalib.cc:449
table::Key< double > calibrationErr
Definition: PhotoCalib.cc:370
table::Key< double > fluxMag0
Definition: PhotoCalib.cc:451
table::Key< double > expTime
Definition: PhotoCalib.cc:450
table::Key< table::Flag > isConstant
Definition: PhotoCalib.cc:371
table::Key< double > fluxMag0Err
Definition: PhotoCalib.cc:452
table::Key< int > version
Definition: PhotoCalib.cc:373
table::Key< int > field
Definition: PhotoCalib.cc:372
Implementation of the Photometric Calibration class.
std::ostream * os
Definition: Schema.cc:557
SchemaMapper * mapper
Definition: SchemaMapper.cc:71
#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:73
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition: MaskedImage.h:1051
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition: MaskedImage.h:1018
The photometric calibration of an exposure.
Definition: PhotoCalib.h:114
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
Definition: BaseRecord.h:80
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:151
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
Key< T > addField(Field< T > const &field, bool doReplace=false)
Add a new field to the Schema, and return the associated Key.
Definition: Schema.cc:479
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.
Definition: SchemaMapper.h:21
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:570
An object passed to Persistable::write to allow it to persist itself.
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.
Definition: Persistable.cc:18
Interface supporting iteration over heterogenous containers.
Definition: Storable.h:58
Class for storing generic metadata.
Definition: PropertySet.h:66
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).
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 find(T... args)
bool strip
Definition: fits.cc:911
T hypot(T... args)
T log(T... args)
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:108
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
std::shared_ptr< PhotoCalib > makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err)
Construct a PhotoCalib from the deprecated Calib-style instFluxMag0/instFluxMag0Err values.
Definition: PhotoCalib.cc:613
bool operator==(FilterProperty const &rhs) const noexcept
Return true iff two FilterProperties are identical.
std::string getPersistenceName() const override
std::shared_ptr< PhotoCalib > makePhotoCalibFromMetadata(daf::base::PropertySet &metadata, bool strip=false)
Construct a PhotoCalib from FITS FLUXMAG0/FLUXMAG0ERR keywords.
Definition: PhotoCalib.cc:595
void write(OutputArchiveHandle &handle) const override
FilterProperty & operator=(FilterProperty const &)=default
std::ostream & operator<<(std::ostream &os, PhotoCalib const &photoCalib)
Definition: PhotoCalib.cc:256
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
Angle abs(Angle const &a)
Definition: Angle.h:106
A base class for image defects.
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
Key< int > photoCalib
Definition: Exposure.cc:67
T to_string(T... args)
Utilities for converting between flux and magnitude in C++.
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override
Definition: warpExposure.cc:0