LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
wcsUtils.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2017 LSST Corporation.
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 <iostream>
24 #include <sstream>
25 #include <cmath>
26 #include <cstring>
27 #include <limits>
28 
29 #include "lsst/daf/base.h"
30 #include "lsst/pex/exceptions.h"
32 #include "lsst/afw/geom/wcsUtils.h"
33 
34 namespace lsst {
35 namespace afw {
36 namespace geom {
37 namespace {
38 /*
39  * Get the name of a SIP matrix coefficient card
40  *
41  * This works for matrix coefficients of the form: <name>_i_j where i and j are zero-based.
42  * Thus it works for SIP matrix terms but not CD matrix terms (because CD terms use 1-based indexing
43  * and have no underscore between the name and the first index).
44  */
45 std::string getSipCoeffCardName(std::string const& name, int i, int j) {
46  return name + std::to_string(i) + "_" + std::to_string(j);
47 }
48 
49 } // namespace
50 
52  lsst::geom::Point2I const& xy0) {
54 
55  wcsMetaData->set("CTYPE1" + wcsName, "LINEAR", "Type of projection");
56  wcsMetaData->set("CTYPE2" + wcsName, "LINEAR", "Type of projection");
57  wcsMetaData->set("CRPIX1" + wcsName, static_cast<double>(1), "Column Pixel Coordinate of Reference");
58  wcsMetaData->set("CRPIX2" + wcsName, static_cast<double>(1), "Row Pixel Coordinate of Reference");
59  wcsMetaData->set("CRVAL1" + wcsName, static_cast<double>(xy0[0]), "Column pixel of Reference Pixel");
60  wcsMetaData->set("CRVAL2" + wcsName, static_cast<double>(xy0[1]), "Row pixel of Reference Pixel");
61  wcsMetaData->set("CUNIT1" + wcsName, "PIXEL", "Column unit");
62  wcsMetaData->set("CUNIT2" + wcsName, "PIXEL", "Row unit");
63 
64  return wcsMetaData;
65 }
66 
68  std::vector<std::string> const names = {"CRPIX1", "CRPIX2", "CRVAL1", "CRVAL2", "CTYPE1",
69  "CTYPE2", "CUNIT1", "CUNIT2", "CD1_1", "CD1_2",
70  "CD2_1", "CD2_2", "WCSAXES"};
71  for (auto const& name : names) {
72  if (metadata.exists(name + wcsName)) {
73  metadata.remove(name + wcsName);
74  }
75  }
76 }
77 
78 Eigen::Matrix2d getCdMatrixFromMetadata(daf::base::PropertySet& metadata) {
79  Eigen::Matrix2d matrix;
80  bool found{false};
81  for (int i = 0; i < 2; ++i) {
82  for (int j = 0; j < 2; ++j) {
83  std::string const cardName = "CD" + std::to_string(i + 1) + "_" + std::to_string(j + 1);
84  if (metadata.exists(cardName)) {
85  matrix(i, j) = metadata.getAsDouble(cardName);
86  found = true;
87  } else {
88  matrix(i, j) = 0.0;
89  }
90  }
91  }
92  if (!found) {
93  throw LSST_EXCEPT(pex::exceptions::TypeError, "No CD matrix coefficients found");
94  }
95  return matrix;
96 }
97 
99  bool strip) {
100  int x0 = 0; // Our value of X0
101  int y0 = 0; // Our value of Y0
102 
103  if (metadata.exists("CRPIX1" + wcsName) && metadata.exists("CRPIX2" + wcsName) &&
104  metadata.exists("CRVAL1" + wcsName) && metadata.exists("CRVAL2" + wcsName) &&
105  (metadata.getAsDouble("CRPIX1" + wcsName) == 1.0) &&
106  (metadata.getAsDouble("CRPIX2" + wcsName) == 1.0)) {
107  x0 = static_cast<int>(std::lround(metadata.getAsDouble("CRVAL1" + wcsName)));
108  y0 = static_cast<int>(std::lround(metadata.getAsDouble("CRVAL2" + wcsName)));
109  if (strip) {
110  deleteBasicWcsMetadata(metadata, wcsName);
111  }
112  }
113  return lsst::geom::Point2I(x0, y0);
114 }
115 
116 Eigen::MatrixXd getSipMatrixFromMetadata(daf::base::PropertySet const& metadata, std::string const& name) {
117  std::string cardName = name + "_ORDER";
118  if (!metadata.exists(cardName)) {
120  "metadata does not contain SIP matrix " + name + ": " + cardName + " not found");
121  };
122  int order = metadata.getAsInt(cardName);
123  if (order < 0) {
125  "matrix order = " + std::to_string(order) + " < 0");
126  }
127  Eigen::MatrixXd matrix(order + 1, order + 1);
128  auto const coeffName = name + "_";
129  for (int i = 0; i <= order; ++i) {
130  for (int j = 0; j <= order; ++j) {
131  std::string const cardName = getSipCoeffCardName(coeffName, i, j);
132  if (metadata.exists(cardName)) {
133  matrix(i, j) = metadata.getAsDouble(cardName);
134  } else {
135  matrix(i, j) = 0.0;
136  }
137  }
138  }
139  return matrix;
140 }
141 
142 bool hasSipMatrix(daf::base::PropertySet const& metadata, std::string const& name) {
143  std::string cardName = name + "_ORDER";
144  if (!metadata.exists(cardName)) {
145  return false;
146  }
147  return metadata.getAsInt(cardName) >= 0;
148 }
149 
151  std::string const& name) {
152  if (matrix.rows() != matrix.cols() || matrix.rows() < 1) {
154  "Matrix must be square and at least 1 x 1; dimensions = " +
155  std::to_string(matrix.rows()) + " x " + std::to_string(matrix.cols()));
156  }
157  int const order = matrix.rows() - 1;
158  auto metadata = std::make_shared<daf::base::PropertyList>();
159  std::string cardName = name + "_ORDER";
160  metadata->set(cardName, order);
161  auto const coeffName = name + "_";
162  for (int i = 0; i <= order; ++i) {
163  for (int j = 0; j <= order; ++j) {
164  if (matrix(i, j) != 0.0) {
165  metadata->set(getSipCoeffCardName(coeffName, i, j), matrix(i, j));
166  }
167  }
168  }
169  return metadata;
170 }
171 
174  Eigen::Matrix2d const& cdMatrix,
175  std::string const& projection) {
176  auto pl = std::make_shared<daf::base::PropertyList>();
177  pl->add("RADESYS", "ICRS");
178  pl->add("CTYPE1", "RA---" + projection);
179  pl->add("CTYPE2", "DEC--" + projection);
180  pl->add("CRPIX1", crpix[0] + 1);
181  pl->add("CRPIX2", crpix[1] + 1);
182  pl->add("CRVAL1", crval[0].asDegrees());
183  pl->add("CRVAL2", crval[1].asDegrees());
184  pl->add("CUNIT1", "deg");
185  pl->add("CUNIT2", "deg");
186  for (int i = 0; i < 2; ++i) {
187  for (int j = 0; j < 2; ++j) {
188  if (cdMatrix(i, j) != 0.0) {
189  std::string name = "CD" + std::to_string(i + 1) + "_" + std::to_string(j + 1);
190  pl->add(name, cdMatrix(i, j));
191  }
192  }
193  }
194  return pl;
195 }
196 
199  Eigen::Matrix2d const& cdMatrix,
200  Eigen::MatrixXd const& sipA,
201  Eigen::MatrixXd const& sipB) {
202  auto metadata = makeSimpleWcsMetadata(crpix, crval, cdMatrix, "TAN-SIP");
203  metadata->combine(makeSipMatrixMetadata(sipA, "A"));
204  metadata->combine(makeSipMatrixMetadata(sipB, "B"));
205  return metadata;
206 }
207 
210  Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA, Eigen::MatrixXd const& sipB,
211  Eigen::MatrixXd const& sipAp, Eigen::MatrixXd const& sipBp) {
212  auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB);
213  metadata->combine(makeSipMatrixMetadata(sipAp, "AP"));
214  metadata->combine(makeSipMatrixMetadata(sipBp, "BP"));
215  return metadata;
216 }
217 
218 } // namespace geom
219 } // namespace afw
220 } // 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::PointKey< double > crpix
Definition: OldWcs.cc:131
table::PointKey< double > crval
Definition: OldWcs.cc:130
Class for storing ordered metadata with comments.
Definition: PropertyList.h:68
Class for storing generic metadata.
Definition: PropertySet.h:67
int getAsInt(std::string const &name) const
Get the last value for a bool/char/short/int property name (possibly hierarchical).
virtual void remove(std::string const &name)
Remove all values for a property name (possibly hierarchical).
bool exists(std::string const &name) const
Determine if a name (possibly hierarchical) exists.
double getAsDouble(std::string const &name) const
Get the last value for any arithmetic property name (possibly hierarchical).
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
Reports errors from accepting an object of an unexpected or inappropriate type.
Definition: Runtime.h:167
bool strip
Definition: fits.cc:911
std::shared_ptr< daf::base::PropertyList > makeSimpleWcsMetadata(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, std::string const &projection="TAN")
Make FITS metadata for a simple FITS WCS (one with no distortion).
Definition: wcsUtils.cc:172
Eigen::MatrixXd getSipMatrixFromMetadata(daf::base::PropertySet const &metadata, std::string const &name)
Definition: wcsUtils.cc:116
bool hasSipMatrix(daf::base::PropertySet const &metadata, std::string const &name)
Definition: wcsUtils.cc:142
std::shared_ptr< daf::base::PropertyList > makeTanSipMetadata(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, Eigen::MatrixXd const &sipA, Eigen::MatrixXd const &sipB)
Make metadata for a TAN-SIP WCS without inverse matrices.
Definition: wcsUtils.cc:197
Eigen::Matrix2d getCdMatrixFromMetadata(daf::base::PropertySet &metadata)
Read a CD matrix from FITS WCS metadata.
Definition: wcsUtils.cc:78
std::shared_ptr< daf::base::PropertyList > createTrivialWcsMetadata(std::string const &wcsName, lsst::geom::Point2I const &xy0)
Definition: wcsUtils.cc:51
std::shared_ptr< daf::base::PropertyList > makeSipMatrixMetadata(Eigen::MatrixXd const &matrix, std::string const &name)
Definition: wcsUtils.cc:150
lsst::geom::Point2I getImageXY0FromMetadata(daf::base::PropertySet &metadata, std::string const &wcsName, bool strip=false)
Definition: wcsUtils.cc:98
void deleteBasicWcsMetadata(daf::base::PropertySet &metadata, std::string const &wcsName)
Definition: wcsUtils.cc:67
Point< int, 2 > Point2I
Definition: Point.h:321
A base class for image defects.
T lround(T... args)
T to_string(T... args)