LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
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 {
41  geom::Box2D bbox;
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  mainSchema.getCitizen().markPersistent();
201  inputDataSchema.getCitizen().markPersistent();
202  }
203 };
204 
205 void CoaddTransmissionCurve::write(OutputArchiveHandle& handle) const {
206  auto const& keys = PersistenceHelper::get();
207  auto mainCat = handle.makeCatalog(keys.mainSchema);
208  auto mainRecord = mainCat.addNew();
209  mainRecord->set(keys.coaddWcs, handle.put(_coaddWcs));
210  mainRecord->set(keys.wavelengthMin, _wavelengthBounds.first);
211  mainRecord->set(keys.wavelengthMax, _wavelengthBounds.second);
212  mainRecord->set(keys.throughputAtMin, _throughputAtBounds.first);
213  mainRecord->set(keys.throughputAtMax, _throughputAtBounds.second);
214  handle.saveCatalog(mainCat);
215  auto inputDataCat = handle.makeCatalog(keys.inputDataSchema);
216  for (auto const& input : _inputs) {
217  auto inputDataRecord = inputDataCat.addNew();
218  inputDataRecord->set(keys.transmission, handle.put(input.transmission));
219  inputDataRecord->set(keys.sensorWcs, handle.put(input.sensorWcs));
220  inputDataRecord->set(keys.validPolygon, handle.put(input.validPolygon));
221  inputDataRecord->set(keys.bbox, input.bbox);
222  inputDataRecord->set(keys.weight, input.weight);
223  }
224  handle.saveCatalog(inputDataCat);
225 }
226 
227 class CoaddTransmissionCurve::Factory : public afw::table::io::PersistableFactory {
228 public:
229  std::shared_ptr<afw::table::io::Persistable> read(InputArchive const& archive,
230  CatalogVector const& catalogs) const override {
231  auto const& keys = PersistenceHelper::get();
232  LSST_ARCHIVE_ASSERT(catalogs.size() == 2u);
233  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.mainSchema);
234  LSST_ARCHIVE_ASSERT(catalogs.back().getSchema() == keys.inputDataSchema);
235  auto const& mainRecord = catalogs.front().front();
237  inputs.reserve(catalogs.back().size());
238  for (auto const& inputDataRecord : catalogs.back()) {
239  CoaddInputData input = {
240  archive.get<afw::image::TransmissionCurve>(inputDataRecord.get(keys.transmission)),
241  archive.get<afw::geom::SkyWcs>(inputDataRecord.get(keys.sensorWcs)),
242  archive.get<afw::geom::polygon::Polygon>(inputDataRecord.get(keys.validPolygon)),
243  inputDataRecord.get(keys.bbox), inputDataRecord.get(keys.weight)};
244  inputs.push_back(std::move(input));
245  }
246  return std::make_shared<CoaddTransmissionCurve>(
247  archive.get<afw::geom::SkyWcs>(mainRecord.get(keys.coaddWcs)),
248  std::make_pair(mainRecord.get(keys.wavelengthMin), mainRecord.get(keys.wavelengthMax)),
249  std::make_pair(mainRecord.get(keys.throughputAtMin), mainRecord.get(keys.throughputAtMax)),
250  std::move(inputs));
251  }
252 
253  Factory(std::string const& name) : afw::table::io::PersistableFactory(name) {}
254 };
255 
256 CoaddTransmissionCurve::Factory CoaddTransmissionCurve::registration("CoaddTransmissionCurve");
257 
258 } // namespace
259 
262  return std::make_shared<CoaddTransmissionCurve>(coaddWcs, inputSensors);
263 }
264 
265 } // namespace algorithms
266 } // namespace meas
267 } // namespace lsst
double weight
BoxKey< lsst::geom::Box2D > Box2DKey
Definition: aggregates.h:202
std::shared_ptr< afw::geom::SkyWcs const > sensorWcs
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.
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
Point< double, 2 > Point2D
Definition: Point.h:324
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
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:168
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.
#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:67
ExposureCatalogT< ExposureRecord > ExposureCatalog
Definition: Exposure.h:499
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