LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
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 
33 #include "lsst/afw/geom/wcsUtils.h"
34 #include "lsst/afw/geom/SkyWcs.h"
35 
36 namespace lsst {
37 namespace afw {
38 namespace geom {
39 namespace {
40 
41 /*
42  * Does one string end with another?
43  *
44  * From https://stackoverflow.com/a/2072890
45  */
46 inline 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  */
60 std::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  */
69 std::shared_ptr<daf::base::PropertyList> getOldSipMetadata(table::BaseRecord const& record);
70 
71 // Unpersist a SkyWcs from table data for an afw::image::Wcs
72 class OldWcsFactory : public table::io::PersistableFactory {
73 public:
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 
86 OldWcsFactory registerWcs("Wcs");
87 
88 // Unpersist a SkyWcs from table data for an afw::image::TanWcs, which may have a record for SIP terms
89 class OldTanWcsFactory : public table::io::PersistableFactory {
90 public:
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 
122 OldTanWcsFactory 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.
126 struct 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 
151 private:
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 
166 std::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 
193 std::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
FilterProperty & operator=(FilterProperty const &)=default
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