LSST Applications g042eb84c57+730a74494b,g04e9c324dd+8c5ae1fdc5,g134cb467dc+1f1e3e7524,g199a45376c+0ba108daf9,g1fd858c14a+fa7d31856b,g210f2d0738+f66ac109ec,g262e1987ae+83a3acc0e5,g29ae962dfc+d856a2cb1f,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+a1e0c9f713,g47891489e3+0d594cb711,g4d44eb3520+c57ec8f3ed,g4d7b6aa1c5+f66ac109ec,g53246c7159+8c5ae1fdc5,g56a1a4eaf3+fd7ad03fde,g64539dfbff+f66ac109ec,g67b6fd64d1+0d594cb711,g67fd3c3899+f66ac109ec,g6985122a63+0d594cb711,g74acd417e5+3098891321,g786e29fd12+668abc6043,g81db2e9a8d+98e2ab9f28,g87389fa792+8856018cbb,g89139ef638+0d594cb711,g8d7436a09f+80fda9ce03,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+a8a29bda4b,g99822b682c+e3ec3c61f9,g9d5c6a246b+0d5dac0c3d,ga41d0fce20+9243b26dd2,gbf99507273+8c5ae1fdc5,gd7ef33dd92+0d594cb711,gdab6d2f7ff+3098891321,ge410e46f29+0d594cb711,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.38
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;
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
170 lsst::geom::SpherePoint const& crval,
171 Eigen::Matrix2d const& cdMatrix,
172 std::string const& projection) {
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
195 lsst::geom::SpherePoint const& crval,
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
206 lsst::geom::Point2D const& crpix, lsst::geom::SpherePoint const& crval,
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
Class for storing ordered metadata with comments.
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
T make_shared(T... args)
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< double, 2 > Point2D
Definition Point.h:324
Point< int, 2 > Point2I
Definition Point.h:321
T lround(T... args)
T to_string(T... args)