LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
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
def write(self, patchRef, catalog)
Write the output.
double weight
BoxKey< lsst::geom::Box2D > Box2DKey
Definition: aggregates.h:202
std::shared_ptr< afw::geom::SkyWcs const > sensorWcs
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
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.
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
afw::table::Key< int > coaddWcs
afw::table::Schema mainSchema
afw::table::Key< double > wavelengthMin
afw::table::Key< double > wavelengthMax
afw::table::Schema inputDataSchema
std::shared_ptr< afw::image::TransmissionCurve const > transmission
STL class.
T min(T... args)
T push_back(T... args)
A base class for image defects.
std::shared_ptr< afw::geom::polygon::Polygon const > validPolygon
afw::table::Key< double > throughputAtMin
T make_pair(T... args)
afw::table::Key< double > throughputAtMax
T max(T... args)
T move(T... args)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
STL class.
Custom catalog class for ExposureRecord/Table.
Definition: Exposure.h:67
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
ExposureCatalogT< ExposureRecord > ExposureCatalog
Definition: Exposure.h:489
static BoxKey addFields(Schema &schema, std::string const &name, std::string const &doc, std::string const &unit)
Add _min_x, _min_y, _max_x, _max_y fields to a Schema, and return a BoxKey that points to them...
Definition: aggregates.cc:60
T reserve(T... args)
geom::Box2D bbox