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