LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
CoaddTransmissionCurve.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2018 LSST/AURA.
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 
24 #include "lsst/geom/Box.h"
25 #include "lsst/afw/geom/SkyWcs.h"
30 
31 namespace lsst {
32 namespace meas {
33 namespace algorithms {
34 
35 namespace {
36 
37 struct CoaddInputData {
42  double weight;
43 };
44 
45 class CoaddTransmissionCurve : public afw::image::TransmissionCurve {
46 public:
48  afw::table::ExposureCatalog const& inputSensors)
49  : _coaddWcs(coaddWcs),
50  _wavelengthBounds(std::numeric_limits<double>::infinity(), // min = +inf, max = -inf *for now*
52  _throughputAtBounds(0.0, 0.0) {
53  _inputs.reserve(inputSensors.size());
54  afw::table::Key<double> weightKey = inputSensors.getSchema()["weight"];
55  double weightSum = 0.0;
56  for (auto const& record : inputSensors) {
57  double const weight = record.get(weightKey);
58  weightSum += weight;
59  auto transmission = record.getTransmissionCurve();
60  if (transmission == nullptr) {
61  throw LSST_EXCEPT(
62  pex::exceptions::InvalidParameterError,
63  (boost::format("No TransmissionCurve for input with ID %d") % record.getId()).str());
64  }
65  if (record.getWcs() == nullptr) {
66  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
67  (boost::format("No Wcs for input with ID %d") % record.getId()).str());
68  }
69  // Set wavelength bounds from the overall minimum and maximum wavelengths from all inputs.
70  auto const inputWavelengthBounds = transmission->getWavelengthBounds();
71  _wavelengthBounds.first = std::min(_wavelengthBounds.first, inputWavelengthBounds.first);
72  _wavelengthBounds.second = std::max(_wavelengthBounds.second, inputWavelengthBounds.second);
73  // Set throughput at and outside bounds from weighted average of throughput at and outside bounds.
74  auto const inputThroughputAtBounds = transmission->getThroughputAtBounds();
75  _throughputAtBounds.first += inputThroughputAtBounds.first * weight;
76  _throughputAtBounds.second += inputThroughputAtBounds.second * weight;
77  // Add an element to the vector with all the stuff we need to store for each epoch.
78  CoaddInputData input = {transmission, record.getWcs(), record.getValidPolygon(),
79  geom::Box2D(record.getBBox()), weight};
80  _inputs.push_back(std::move(input));
81  }
82  _throughputAtBounds.first /= weightSum;
83  _throughputAtBounds.second /= weightSum;
84  }
85 
86  // Private constructor used only for persistence
87  CoaddTransmissionCurve(std::shared_ptr<afw::geom::SkyWcs> coaddWcs,
88  std::pair<double, double> wavelengthBounds,
89  std::pair<double, double> throughputAtBounds, std::vector<CoaddInputData> inputs)
90  : _coaddWcs(std::move(coaddWcs)),
91  _wavelengthBounds(std::move(wavelengthBounds)),
92  _throughputAtBounds(std::move(throughputAtBounds)),
93  _inputs(std::move(inputs)) {}
94 
95  // All TransmissionCurves are immutable and noncopyable.
96  CoaddTransmissionCurve(CoaddTransmissionCurve const&) = delete;
97  CoaddTransmissionCurve(CoaddTransmissionCurve&&) = delete;
98  CoaddTransmissionCurve& operator=(CoaddTransmissionCurve const&) = delete;
99  CoaddTransmissionCurve& operator=(CoaddTransmissionCurve&&) = delete;
100 
101  std::pair<double, double> getWavelengthBounds() const override { return _wavelengthBounds; }
102 
103  std::pair<double, double> getThroughputAtBounds() const override { return _throughputAtBounds; }
104 
105  void sampleAt(geom::Point2D const& position, ndarray::Array<double const, 1, 1> const& wavelengths,
106  ndarray::Array<double, 1, 1> const& out) const override {
107  auto const coord = _coaddWcs->pixelToSky(position);
108  ndarray::Array<double, 1, 1> tmp = ndarray::allocate(out.getShape());
109  tmp.deep() = 0.0;
110  out.deep() = 0.0;
111  double weightSum = 0.0;
112  for (auto const& input : _inputs) {
113  geom::Point2D const inputPosition = input.sensorWcs->skyToPixel(coord);
114  if (!input.bbox.contains(inputPosition)) {
115  continue;
116  }
117  if (input.validPolygon && !input.validPolygon->contains(inputPosition)) {
118  continue;
119  }
120  // note that `tmp` is an output argument here
121  input.transmission->sampleAt(inputPosition, wavelengths, tmp);
122  tmp.deep() *= input.weight;
123  out.deep() += tmp;
124  weightSum += input.weight;
125  }
126  if (weightSum == 0.0) {
127  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
128  (boost::format("No input TransmissionCurves at point (%s, %s)") %
129  position.getX() % position.getY())
130  .str());
131  }
132  out.deep() /= weightSum;
133  }
134 
135  bool isPersistable() const noexcept override {
136  for (auto const& input : _inputs) {
137  if (!input.transmission->isPersistable()) {
138  return false;
139  }
140  }
141  return true;
142  }
143 
144 private:
145  std::string getPersistenceName() const override { return "CoaddTransmissionCurve"; }
146 
147  std::string getPythonModule() const override { return "lsst.meas.algorithms"; }
148 
149  struct PersistenceHelper;
150 
151  class Factory;
152 
153  void write(OutputArchiveHandle& handle) const override;
154 
155  static Factory registration;
156 
158  std::pair<double, double> _wavelengthBounds;
159  std::pair<double, double> _throughputAtBounds;
161 };
162 
163 struct CoaddTransmissionCurve::PersistenceHelper {
164  afw::table::Schema mainSchema;
165  afw::table::Key<int> coaddWcs;
166  afw::table::Key<double> wavelengthMin;
167  afw::table::Key<double> wavelengthMax;
168  afw::table::Key<double> throughputAtMin;
169  afw::table::Key<double> throughputAtMax;
170  afw::table::Schema inputDataSchema;
171  afw::table::Key<int> transmission;
172  afw::table::Key<int> sensorWcs;
173  afw::table::Key<int> validPolygon;
175  afw::table::Key<double> weight;
176 
177  static PersistenceHelper const& get() {
178  static PersistenceHelper const instance;
179  return instance;
180  }
181 
182 private:
183  PersistenceHelper()
184  : mainSchema(),
185  coaddWcs(mainSchema.addField<int>("coaddWcs", "archive ID for the coadd's WCS")),
186  wavelengthMin(mainSchema.addField<double>("wavelengthMin", "getWavelengthBounds().min")),
187  wavelengthMax(mainSchema.addField<double>("wavelengthMax", "getWavelengthBounds().max")),
188  throughputAtMin(mainSchema.addField<double>("throughputAtMin", "throughputAtBounds().first")),
189  throughputAtMax(mainSchema.addField<double>("throughputAtMax", "throughputAtBounds().second")),
190  inputDataSchema(),
191  transmission(inputDataSchema.addField<int>("transmission",
192  "archive ID for this sensor's TransmissionCurve")),
193  sensorWcs(inputDataSchema.addField<int>("sensorWcs", "archive ID for this sensor's WCS")),
194  validPolygon(inputDataSchema.addField<int>("validPolygon",
195  "archive ID for this sensor's ValidPolygon")),
196  bbox(afw::table::Box2DKey::addFields(inputDataSchema, "bbox", "bounding box of the sensor",
197  "pixel")),
198  weight(inputDataSchema.addField<double>("weight",
199  "relative weight for this sensor in the average")) {}
200 };
201 
202 void CoaddTransmissionCurve::write(OutputArchiveHandle& handle) const {
203  auto const& keys = PersistenceHelper::get();
204  auto mainCat = handle.makeCatalog(keys.mainSchema);
205  auto mainRecord = mainCat.addNew();
206  mainRecord->set(keys.coaddWcs, handle.put(_coaddWcs));
207  mainRecord->set(keys.wavelengthMin, _wavelengthBounds.first);
208  mainRecord->set(keys.wavelengthMax, _wavelengthBounds.second);
209  mainRecord->set(keys.throughputAtMin, _throughputAtBounds.first);
210  mainRecord->set(keys.throughputAtMax, _throughputAtBounds.second);
211  handle.saveCatalog(mainCat);
212  auto inputDataCat = handle.makeCatalog(keys.inputDataSchema);
213  for (auto const& input : _inputs) {
214  auto inputDataRecord = inputDataCat.addNew();
215  inputDataRecord->set(keys.transmission, handle.put(input.transmission));
216  inputDataRecord->set(keys.sensorWcs, handle.put(input.sensorWcs));
217  inputDataRecord->set(keys.validPolygon, handle.put(input.validPolygon));
218  inputDataRecord->set(keys.bbox, input.bbox);
219  inputDataRecord->set(keys.weight, input.weight);
220  }
221  handle.saveCatalog(inputDataCat);
222 }
223 
224 class CoaddTransmissionCurve::Factory : public afw::table::io::PersistableFactory {
225 public:
226  std::shared_ptr<afw::table::io::Persistable> read(InputArchive const& archive,
227  CatalogVector const& catalogs) const override {
228  auto const& keys = PersistenceHelper::get();
229  LSST_ARCHIVE_ASSERT(catalogs.size() == 2u);
230  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.mainSchema);
231  LSST_ARCHIVE_ASSERT(catalogs.back().getSchema() == keys.inputDataSchema);
232  auto const& mainRecord = catalogs.front().front();
234  inputs.reserve(catalogs.back().size());
235  for (auto const& inputDataRecord : catalogs.back()) {
236  CoaddInputData input = {
237  archive.get<afw::image::TransmissionCurve>(inputDataRecord.get(keys.transmission)),
238  archive.get<afw::geom::SkyWcs>(inputDataRecord.get(keys.sensorWcs)),
239  archive.get<afw::geom::polygon::Polygon>(inputDataRecord.get(keys.validPolygon)),
240  inputDataRecord.get(keys.bbox), inputDataRecord.get(keys.weight)};
241  inputs.push_back(std::move(input));
242  }
243  return std::make_shared<CoaddTransmissionCurve>(
244  archive.get<afw::geom::SkyWcs>(mainRecord.get(keys.coaddWcs)),
245  std::make_pair(mainRecord.get(keys.wavelengthMin), mainRecord.get(keys.wavelengthMax)),
246  std::make_pair(mainRecord.get(keys.throughputAtMin), mainRecord.get(keys.throughputAtMax)),
247  std::move(inputs));
248  }
249 
250  Factory(std::string const& name) : afw::table::io::PersistableFactory(name) {}
251 };
252 
253 CoaddTransmissionCurve::Factory CoaddTransmissionCurve::registration("CoaddTransmissionCurve");
254 
255 } // namespace
256 
259  return std::make_shared<CoaddTransmissionCurve>(coaddWcs, inputSensors);
260 }
261 
262 } // namespace algorithms
263 } // namespace meas
264 } // namespace lsst
table::Key< std::string > name
Definition: Amplifier.cc:116
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
Custom catalog class for ExposureRecord/Table.
Definition: Exposure.h:311
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
T make_pair(T... args)
T max(T... args)
T min(T... args)
T move(T... args)
std::string getPythonModule() const override
std::string getPersistenceName() const override
void write(OutputArchiveHandle &handle) const override
FilterProperty & operator=(FilterProperty const &)=default
BoxKey< lsst::geom::Box2D > Box2DKey
Definition: aggregates.h:202
ExposureCatalogT< ExposureRecord > ExposureCatalog
Definition: Exposure.h:489
std::shared_ptr< afw::image::TransmissionCurve const > makeCoaddTransmissionCurve(std::shared_ptr< afw::geom::SkyWcs const > coaddWcs, afw::table::ExposureCatalog const &inputSensors)
Create a TransmissionCurve that represents the effective throughput on a coadd.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
STL namespace.
T push_back(T... args)
T reserve(T... args)
double weight
afw::table::Key< int > coaddWcs
std::shared_ptr< afw::geom::polygon::Polygon const > validPolygon
std::shared_ptr< afw::geom::SkyWcs const > sensorWcs
std::shared_ptr< afw::image::TransmissionCurve const > transmission
geom::Box2D bbox
afw::table::Key< double > wavelengthMin
afw::table::Key< double > wavelengthMax
afw::table::Schema inputDataSchema
afw::table::Schema mainSchema
afw::table::Key< double > throughputAtMin
afw::table::Key< double > throughputAtMax