LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"
30
31namespace lsst {
32namespace meas {
33namespace algorithms {
34
35namespace {
36
37struct CoaddInputData {
42 double weight;
43};
44
45class CoaddTransmissionCurve : public afw::image::TransmissionCurve {
46public:
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,
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
144private:
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
163struct 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
182private:
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")),
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
202void 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
224class CoaddTransmissionCurve::Factory : public afw::table::io::PersistableFactory {
225public:
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
253CoaddTransmissionCurve::Factory CoaddTransmissionCurve::registration("CoaddTransmissionCurve");
254
255} // namespace
256
261
262} // namespace algorithms
263} // namespace meas
264} // namespace lsst
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
table::Key< double > throughputAtMax
table::Key< double > throughputAtMin
#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)
ExposureCatalogT< ExposureRecord > ExposureCatalog
Definition Exposure.h:489
BoxKey< lsst::geom::Box2D > Box2DKey
Definition aggregates.h:284
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.
STL namespace.
afw::table::Key< double > weight
afw::table::Key< int > coaddWcs
std::shared_ptr< afw::geom::SkyWcs const > sensorWcs
std::shared_ptr< afw::image::TransmissionCurve const > transmission
afw::table::Key< double > wavelengthMin
afw::table::Key< double > wavelengthMax
afw::table::Schema inputDataSchema
afw::table::Schema mainSchema
Key< int > validPolygon
Definition Exposure.cc:69
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override