LSST Applications g063fba187b+fee0456c91,g0f08755f38+ea96e5a5a3,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+90257ff92a,g20f6ffc8e0+ea96e5a5a3,g217e2c1bcf+937a289c59,g28da252d5a+daa7da44eb,g2bbee38e9b+253935c60e,g2bc492864f+253935c60e,g3156d2b45e+6e55a43351,g32e5bea42b+31359a2a7a,g347aa1857d+253935c60e,g35bb328faa+a8ce1bb630,g3a166c0a6a+253935c60e,g3b1af351f3+a8ce1bb630,g3e281a1b8c+c5dd892a6c,g414038480c+416496e02f,g41af890bb2+afe91b1188,g599934f4f4+0db33f7991,g7af13505b9+e36de7bce6,g80478fca09+da231ba887,g82479be7b0+a4516e59e3,g858d7b2824+ea96e5a5a3,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+bc6ab8dfbd,gb58c049af0+d64f4d3760,gc28159a63d+253935c60e,gcab2d0539d+3f2b72788c,gcf0d15dbbd+4ea9c45075,gda6a2b7d83+4ea9c45075,gdaeeff99f8+1711a396fd,ge79ae78c31+253935c60e,gef2f8181fd+3031e3cf99,gf0baf85859+c1f95f4921,gfa517265be+ea96e5a5a3,gfa999e8aa5+17cd334064,w.2024.50
LSST Data Management Base Package
Loading...
Searching...
No Matches
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 <cmath>
25
26#include "lsst/daf/base.h"
27#include "lsst/pex/exceptions.h"
30
31namespace lsst {
32namespace afw {
33namespace geom {
34namespace {
35/*
36 * Get the name of a SIP matrix coefficient card
37 *
38 * This works for matrix coefficients of the form: <name>_i_j where i and j are zero-based.
39 * Thus it works for SIP matrix terms but not CD matrix terms (because CD terms use 1-based indexing
40 * and have no underscore between the name and the first index).
41 */
42std::string getSipCoeffCardName(std::string const& name, int i, int j) {
43 return name + std::to_string(i) + "_" + std::to_string(j);
44}
45
46} // namespace
47
49 lsst::geom::Point2I const& xy0) {
51
52 wcsMetaData->set("CTYPE1" + wcsName, "LINEAR", "Type of projection");
53 wcsMetaData->set("CTYPE2" + wcsName, "LINEAR", "Type of projection");
54 wcsMetaData->set("CRPIX1" + wcsName, static_cast<double>(1), "Column Pixel Coordinate of Reference");
55 wcsMetaData->set("CRPIX2" + wcsName, static_cast<double>(1), "Row Pixel Coordinate of Reference");
56 wcsMetaData->set("CRVAL1" + wcsName, static_cast<double>(xy0[0]), "Column pixel of Reference Pixel");
57 wcsMetaData->set("CRVAL2" + wcsName, static_cast<double>(xy0[1]), "Row pixel of Reference Pixel");
58 wcsMetaData->set("CUNIT1" + wcsName, "PIXEL", "Column unit");
59 wcsMetaData->set("CUNIT2" + wcsName, "PIXEL", "Row unit");
60
61 return wcsMetaData;
62}
63
65 std::vector<std::string> const names = {"CRPIX1", "CRPIX2", "CRVAL1", "CRVAL2", "CTYPE1",
66 "CTYPE2", "CUNIT1", "CUNIT2", "CD1_1", "CD1_2",
67 "CD2_1", "CD2_2", "WCSAXES"};
68 for (auto const& name : names) {
69 if (metadata.exists(name + wcsName)) {
70 metadata.remove(name + wcsName);
71 }
72 }
73}
74
76 Eigen::Matrix2d matrix;
77 bool found{false};
78 for (int i = 0; i < 2; ++i) {
79 for (int j = 0; j < 2; ++j) {
80 std::string const cardName = "CD" + std::to_string(i + 1) + "_" + std::to_string(j + 1);
81 if (metadata.exists(cardName)) {
82 matrix(i, j) = metadata.getAsDouble(cardName);
83 found = true;
84 } else {
85 matrix(i, j) = 0.0;
86 }
87 }
88 }
89 if (!found) {
90 throw LSST_EXCEPT(pex::exceptions::TypeError, "No CD matrix coefficients found");
91 }
92 return matrix;
93}
94
96 bool strip) {
97 int x0 = 0; // Our value of X0
98 int y0 = 0; // Our value of Y0
99
100 if (metadata.exists("CRPIX1" + wcsName) && metadata.exists("CRPIX2" + wcsName) &&
101 metadata.exists("CRVAL1" + wcsName) && metadata.exists("CRVAL2" + wcsName) &&
102 (metadata.getAsDouble("CRPIX1" + wcsName) == 1.0) &&
103 (metadata.getAsDouble("CRPIX2" + wcsName) == 1.0)) {
104 x0 = static_cast<int>(std::lround(metadata.getAsDouble("CRVAL1" + wcsName)));
105 y0 = static_cast<int>(std::lround(metadata.getAsDouble("CRVAL2" + wcsName)));
106 if (strip) {
107 deleteBasicWcsMetadata(metadata, wcsName);
108 }
109 }
110 return lsst::geom::Point2I(x0, y0);
111}
112
113Eigen::MatrixXd getSipMatrixFromMetadata(daf::base::PropertySet const& metadata, std::string const& name) {
114 std::string cardName = name + "_ORDER";
115 if (!metadata.exists(cardName)) {
117 "metadata does not contain SIP matrix " + name + ": " + cardName + " not found");
118 };
119 int order = metadata.getAsInt(cardName);
120 if (order < 0) {
122 "matrix order = " + std::to_string(order) + " < 0");
123 }
124 Eigen::MatrixXd matrix(order + 1, order + 1);
125 auto const coeffName = name + "_";
126 for (int i = 0; i <= order; ++i) {
127 for (int j = 0; j <= order; ++j) {
128 std::string const cardName = getSipCoeffCardName(coeffName, i, j);
129 if (metadata.exists(cardName)) {
130 matrix(i, j) = metadata.getAsDouble(cardName);
131 } else {
132 matrix(i, j) = 0.0;
133 }
134 }
135 }
136 return matrix;
137}
138
139bool hasSipMatrix(daf::base::PropertySet const& metadata, std::string const& name) {
140 std::string cardName = name + "_ORDER";
141 if (!metadata.exists(cardName)) {
142 return false;
143 }
144 return metadata.getAsInt(cardName) >= 0;
145}
146
148 std::string const& name) {
149 if (matrix.rows() != matrix.cols() || matrix.rows() < 1) {
151 "Matrix must be square and at least 1 x 1; dimensions = " +
152 std::to_string(matrix.rows()) + " x " + std::to_string(matrix.cols()));
153 }
154 int const order = matrix.rows() - 1;
155 auto metadata = std::make_shared<daf::base::PropertyList>();
156 std::string cardName = name + "_ORDER";
157 metadata->set(cardName, order);
158 auto const coeffName = name + "_";
159 for (int i = 0; i <= order; ++i) {
160 for (int j = 0; j <= order; ++j) {
161 if (matrix(i, j) != 0.0) {
162 metadata->set(getSipCoeffCardName(coeffName, i, j), matrix(i, j));
163 }
164 }
165 }
166 return metadata;
167}
168
171 Eigen::Matrix2d const& cdMatrix,
172 std::string const& projection) {
173 auto pl = std::make_shared<daf::base::PropertyList>();
174 pl->add("RADESYS", "ICRS");
175 pl->add("CTYPE1", "RA---" + projection);
176 pl->add("CTYPE2", "DEC--" + projection);
177 pl->add("CRPIX1", crpix[0] + 1);
178 pl->add("CRPIX2", crpix[1] + 1);
179 pl->add("CRVAL1", crval[0].asDegrees());
180 pl->add("CRVAL2", crval[1].asDegrees());
181 pl->add("CUNIT1", "deg");
182 pl->add("CUNIT2", "deg");
183 for (int i = 0; i < 2; ++i) {
184 for (int j = 0; j < 2; ++j) {
185 if (cdMatrix(i, j) != 0.0) {
186 std::string name = "CD" + std::to_string(i + 1) + "_" + std::to_string(j + 1);
187 pl->add(name, cdMatrix(i, j));
188 }
189 }
190 }
191 return pl;
192}
193
196 Eigen::Matrix2d const& cdMatrix,
197 Eigen::MatrixXd const& sipA,
198 Eigen::MatrixXd const& sipB) {
199 auto metadata = makeSimpleWcsMetadata(crpix, crval, cdMatrix, "TAN-SIP");
200 metadata->combine(*makeSipMatrixMetadata(sipA, "A"));
201 metadata->combine(*makeSipMatrixMetadata(sipB, "B"));
202 return metadata;
203}
204
207 Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA, Eigen::MatrixXd const& sipB,
208 Eigen::MatrixXd const& sipAp, Eigen::MatrixXd const& sipBp) {
209 auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB);
210 metadata->combine(*makeSipMatrixMetadata(sipAp, "AP"));
211 metadata->combine(*makeSipMatrixMetadata(sipBp, "BP"));
212 return metadata;
213}
214
218
219} // namespace geom
220} // namespace afw
221} // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
table::PointKey< double > crpix
Definition OldWcs.cc:129
table::PointKey< double > crval
Definition OldWcs.cc:128
Class for storing ordered metadata with comments.
Class for storing generic metadata.
Definition PropertySet.h:66
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:930
void stripWcsMetadata(daf::base::PropertySet &metadata)
Strip all WCS metadata from a header.
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:169
Eigen::MatrixXd getSipMatrixFromMetadata(daf::base::PropertySet const &metadata, std::string const &name)
Definition wcsUtils.cc:113
bool hasSipMatrix(daf::base::PropertySet const &metadata, std::string const &name)
Definition wcsUtils.cc:139
void stripWcsMetadata(daf::base::PropertySet &metadata)
Definition wcsUtils.cc:215
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:194
Eigen::Matrix2d getCdMatrixFromMetadata(daf::base::PropertySet &metadata)
Read a CD matrix from FITS WCS metadata.
Definition wcsUtils.cc:75
std::shared_ptr< daf::base::PropertyList > createTrivialWcsMetadata(std::string const &wcsName, lsst::geom::Point2I const &xy0)
Definition wcsUtils.cc:48
std::shared_ptr< daf::base::PropertyList > makeSipMatrixMetadata(Eigen::MatrixXd const &matrix, std::string const &name)
Definition wcsUtils.cc:147
lsst::geom::Point2I getImageXY0FromMetadata(daf::base::PropertySet &metadata, std::string const &wcsName, bool strip=false)
Definition wcsUtils.cc:95
void deleteBasicWcsMetadata(daf::base::PropertySet &metadata, std::string const &wcsName)
Definition wcsUtils.cc:64
Point< int, 2 > Point2I
Definition Point.h:321
T lround(T... args)
T to_string(T... args)
table::Key< int > order