LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
OldWcs.cc
Go to the documentation of this file.
1/*
2 * LSST Data Management System
3 * Copyright 2017 AURA/LSST.
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 <memory>
24
25#include "astshim.h"
26
35
36namespace lsst {
37namespace afw {
38namespace geom {
39namespace {
40
41/*
42 * Does one string end with another?
43 *
44 * From https://stackoverflow.com/a/2072890
45 */
46inline bool endsWith(std::string const& value, std::string const& suffix) {
47 if (suffix.size() > value.size()) {
48 return false;
49 }
50 return std::equal(suffix.rbegin(), suffix.rend(), value.rbegin());
51}
52
53/*
54 * @internal Get FITS WCS metadata from a record for an lsst::afw::image::Wcs
55 * or the non-SIP terms of an lsst::afw::image::TanWcs
56 *
57 * @param[in] record Record holding data for an old lsst::afw::image::Wcs
58 * @return FITS WCS metadata
59 */
60std::shared_ptr<daf::base::PropertyList> getOldWcsMetadata(table::BaseRecord const& record);
61
62/*
63 * @internal Get FITS WCS metadata for the SIP terms of a TAN-SIP WCS
64 * saved by lsst::afw::image::TanWcs
65 *
66 * @param[in] record Record holding SIP terms for an old lsst::afw::image::TanWcs
67 * @return FITS WCS metadata
68 */
69std::shared_ptr<daf::base::PropertyList> getOldSipMetadata(table::BaseRecord const& record);
70
71// Unpersist a SkyWcs from table data for an afw::image::Wcs
72class OldWcsFactory : public table::io::PersistableFactory {
73public:
74 explicit OldWcsFactory(std::string const& name) : table::io::PersistableFactory(name) {}
75
76 std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
77 CatalogVector const& catalogs) const override {
78 LSST_ARCHIVE_ASSERT(catalogs.size() >= 1u);
79 LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
80 auto const& record = catalogs.front().front();
81 auto const metadata = getOldWcsMetadata(record);
82 return std::make_shared<SkyWcs>(*metadata);
83 }
84};
85
86OldWcsFactory registerWcs("Wcs");
87
88// Unpersist a SkyWcs from table data for an afw::image::TanWcs, which may have a record for SIP terms
89class OldTanWcsFactory : public table::io::PersistableFactory {
90public:
91 explicit OldTanWcsFactory(std::string const& name) : table::io::PersistableFactory(name) {}
92
93 std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
94 CatalogVector const& catalogs) const override {
95 LSST_ARCHIVE_ASSERT(catalogs.size() >= 1u);
96 auto const& record = catalogs.front().front();
97 auto const metadata = getOldWcsMetadata(record);
98
99 if (catalogs.size() > 1u) {
100 LSST_ARCHIVE_ASSERT(catalogs.size() == 2u);
101 LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
102 LSST_ARCHIVE_ASSERT(catalogs.back().size() == 1u);
103 std::shared_ptr<table::BaseRecord const> sipRecord = catalogs.back().begin();
104
105 // if CTYPE1 or CTYPE2 does not end with -SIP (and likely it never does), append it
106 std::string const ctypeSuffix = "-SIP";
107 for (auto i = 1; i <= 2; ++i) {
108 std::string const ctypeName = "CTYPE" + std::to_string(i);
109 auto const ctypeValue = metadata->getAsString(ctypeName);
110 if (!endsWith(ctypeValue, ctypeSuffix)) {
111 metadata->set(ctypeName, ctypeValue + ctypeSuffix);
112 }
113 }
114
115 auto sipMetadata = getOldSipMetadata(*sipRecord);
116 metadata->combine(*sipMetadata);
117 }
118 return std::make_shared<SkyWcs>(*metadata);
119 }
120};
121
122OldTanWcsFactory registerTanWcs("TanWcs");
123
124// Read-only singleton struct containing the schema and keys that lsst::afw::image::Wcs
125// and the non-SIP portion of lsst::afw::image::TanWcs were mapped to.
126struct OldWcsPersistenceHelper {
127 table::Schema schema;
128 table::PointKey<double> crval;
129 table::PointKey<double> crpix;
130 table::Key<table::Array<double>> cd;
131 table::Key<std::string> ctype1;
132 table::Key<std::string> ctype2;
133 table::Key<double> equinox;
134 table::Key<std::string> radesys;
135 table::Key<std::string> cunit1;
136 table::Key<std::string> cunit2;
137
138 static OldWcsPersistenceHelper const& get() {
139 static OldWcsPersistenceHelper instance;
140 return instance;
141 };
142
143 // No copying
144 OldWcsPersistenceHelper(const OldWcsPersistenceHelper&) = delete;
145 OldWcsPersistenceHelper& operator=(const OldWcsPersistenceHelper&) = delete;
146
147 // No moving
148 OldWcsPersistenceHelper(OldWcsPersistenceHelper&&) = delete;
149 OldWcsPersistenceHelper& operator=(OldWcsPersistenceHelper&&) = delete;
150
151private:
152 OldWcsPersistenceHelper()
153 : schema(),
154 crval(table::PointKey<double>::addFields(schema, "crval", "celestial reference point", "deg")),
155 crpix(table::PointKey<double>::addFields(schema, "crpix", "pixel reference point", "pixel")),
156 cd(schema.addField<table::Array<double>>(
157 "cd", "linear transform matrix, ordered (1_1, 2_1, 1_2, 2_2)", 4)),
158 ctype1(schema.addField<std::string>("ctype1", "coordinate type", 72)),
159 ctype2(schema.addField<std::string>("ctype2", "coordinate type", 72)),
160 equinox(schema.addField<double>("equinox", "equinox of coordinates")),
161 radesys(schema.addField<std::string>("radesys", "coordinate system for equinox", 72)),
162 cunit1(schema.addField<std::string>("cunit1", "coordinate units", 72)),
163 cunit2(schema.addField<std::string>("cunit2", "coordinate units", 72)) {}
164};
165
166std::shared_ptr<daf::base::PropertyList> getOldWcsMetadata(table::BaseRecord const& record) {
167 auto const& keys = OldWcsPersistenceHelper::get();
168 LSST_ARCHIVE_ASSERT(record.getSchema() == keys.schema);
169
170 auto metadata = std::make_shared<daf::base::PropertyList>();
171
172 auto crvalDeg = record.get(keys.crval);
173 auto crpix = record.get(keys.crpix);
174 auto cd = record.get(keys.cd);
175 metadata->set("CRVAL1", crvalDeg[0]);
176 metadata->set("CRVAL2", crvalDeg[1]);
177 // Add 1 to CRPIX because the field was saved using the LSST standard: 0,0 is lower left pixel
178 metadata->set("CRPIX1", crpix[0] + 1);
179 metadata->set("CRPIX2", crpix[1] + 1);
180 metadata->set("CD1_1", cd[0]);
181 metadata->set("CD2_1", cd[1]);
182 metadata->set("CD1_2", cd[2]);
183 metadata->set("CD2_2", cd[3]);
184 metadata->set("CTYPE1", record.get(keys.ctype1));
185 metadata->set("CTYPE2", record.get(keys.ctype2));
186 metadata->set("EQUINOX", record.get(keys.equinox));
187 metadata->set("RADESYS", record.get(keys.radesys));
188 metadata->set("CUNIT1", record.get(keys.cunit1));
189 metadata->set("CUNIT2", record.get(keys.cunit2));
190 return metadata;
191}
192
193std::shared_ptr<daf::base::PropertyList> getOldSipMetadata(table::BaseRecord const& record) {
194 // Cannot use a PersistenceHelper for the SIP terms because the length of each SIP array
195 // must be read from the schema
196 afw::table::Key<table::Array<double>> kA;
197 afw::table::Key<table::Array<double>> kB;
198 afw::table::Key<table::Array<double>> kAp;
199 afw::table::Key<table::Array<double>> kBp;
200 try {
201 kA = record.getSchema()["A"];
202 kB = record.getSchema()["B"];
203 kAp = record.getSchema()["Ap"];
204 kBp = record.getSchema()["Bp"];
205 } catch (...) {
206 throw LSST_EXCEPT(afw::table::io::MalformedArchiveError,
207 "Incorrect schema for TanWcs distortion terms");
208 }
209
210 // Adding 0.5 and truncating the result guarantees we'll get the right answer
211 // for small ints even when round-off error is involved.
212 int nA = static_cast<int>(std::sqrt(kA.getSize() + 0.5));
213 int nB = static_cast<int>(std::sqrt(kB.getSize() + 0.5));
214 int nAp = static_cast<int>(std::sqrt(kAp.getSize() + 0.5));
215 int nBp = static_cast<int>(std::sqrt(kBp.getSize() + 0.5));
216 if (nA * nA != kA.getSize()) {
217 throw LSST_EXCEPT(table::io::MalformedArchiveError, "Forward X SIP matrix is not square.");
218 }
219 if (nB * nB != kB.getSize()) {
220 throw LSST_EXCEPT(table::io::MalformedArchiveError, "Forward Y SIP matrix is not square.");
221 }
222 if (nAp * nAp != kAp.getSize()) {
223 throw LSST_EXCEPT(table::io::MalformedArchiveError, "Reverse X SIP matrix is not square.");
224 }
225 if (nBp * nBp != kBp.getSize()) {
226 throw LSST_EXCEPT(table::io::MalformedArchiveError, "Reverse Y SIP matrix is not square.");
227 }
228 Eigen::Map<Eigen::MatrixXd const> mapA((record)[kA].getData(), nA, nA);
229 Eigen::Map<Eigen::MatrixXd const> mapB((record)[kB].getData(), nB, nB);
230 Eigen::Map<Eigen::MatrixXd const> mapAp((record)[kAp].getData(), nAp, nAp);
231 Eigen::Map<Eigen::MatrixXd const> mapBp((record)[kBp].getData(), nBp, nBp);
232
233 auto metadata = makeSipMatrixMetadata(mapA, "A");
234 metadata->combine(*makeSipMatrixMetadata(mapB, "B"));
235 metadata->combine(*makeSipMatrixMetadata(mapAp, "AP"));
236 metadata->combine(*makeSipMatrixMetadata(mapBp, "BP"));
237 return metadata;
238}
239
240} // namespace
241} // namespace geom
242} // namespace afw
243} // 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
table::Key< std::string > ctype1
Definition: OldWcs.cc:131
table::PointKey< double > crpix
Definition: OldWcs.cc:129
table::Key< std::string > radesys
Definition: OldWcs.cc:134
table::Schema schema
Definition: OldWcs.cc:127
table::Key< table::Array< double > > cd
Definition: OldWcs.cc:130
table::Key< std::string > cunit2
Definition: OldWcs.cc:136
table::Key< std::string > cunit1
Definition: OldWcs.cc:135
table::Key< std::string > ctype2
Definition: OldWcs.cc:132
table::Key< double > equinox
Definition: OldWcs.cc:133
table::PointKey< double > crval
Definition: OldWcs.cc:128
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
T equal(T... args)
std::shared_ptr< daf::base::PropertyList > makeSipMatrixMetadata(Eigen::MatrixXd const &matrix, std::string const &name)
Definition: wcsUtils.cc:147
A base class for image defects.
STL namespace.
T rbegin(T... args)
T rend(T... args)
T size(T... args)
T sqrt(T... args)
T to_string(T... args)
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override
Definition: warpExposure.cc:0