LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
SpherePoint.cc
Go to the documentation of this file.
1/*
2 * Developed for the LSST Data Management System.
3 * This product includes software developed by the LSST Project
4 * (https://www.lsst.org).
5 * See the COPYRIGHT file at the top-level directory of this distribution
6 * for details of code ownership.
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 GNU General Public License
19 * along with this program. If not, see <https://www.gnu.org/licenses/>.
20 */
21
22#include <cmath>
23#include <string>
24#include <sstream>
25
26#include "Eigen/Geometry"
27
31#include "lsst/pex/exceptions.h"
32
34using namespace std;
35
36namespace lsst {
37
38namespace geom {
40namespace {
54Angle haversine(Angle const& deltaLon, Angle const& deltaLat, double cosLat1, double cosLat2) {
55 double const sinLat = sin(deltaLat / 2.0);
56 double const sinLon = sin(deltaLon / 2.0);
57 double const havDDelta = sinLat * sinLat;
58 double const havDAlpha = sinLon * sinLon;
59 double const havD = havDDelta + cosLat1 * cosLat2 * havDAlpha;
60 double const sinDHalf = sqrt(havD);
61 return (2.0 * asin(sinDHalf)) * radians;
62}
63} // end namespace
64
66 : SpherePoint(longitude * units, latitude * units) {}
67
69 : _longitude(longitude.wrap().asRadians()), _latitude(latitude.asRadians()) {
70 if (fabs(_latitude) > HALFPI) {
71 throw pexExcept::InvalidParameterError("Angle " + to_string(latitude.asDegrees()) +
72 " is not a valid latitude.");
73 }
74}
75
77 // sphgeom Vector3d has its own normalization,
78 // but its behavior is not documented for non-finite values
79 double norm = vector.getNorm();
80 if (norm <= 0.0) {
81 stringstream buffer;
82 buffer << "Vector " << vector << " has zero norm and cannot be normalized.";
84 }
85 // To avoid unexpected behavior from mixing finite elements and infinite norm
86 if (!isfinite(norm)) {
87 norm = NAN;
88 }
89 auto unitVector =
91 _set(unitVector);
92}
93
95 : SpherePoint(lonLat.getLon().asRadians(), lonLat.getLat().asRadians(), radians) {}
96
97void SpherePoint::_set(sphgeom::UnitVector3d const& unitVector) {
98 _latitude = asin(unitVector.z());
99 if (!atPole()) {
100 // Need to convert to Angle, Angle::wrap, and convert back to radians
101 // to handle _longitude = -1e-16 without code duplication
102 _longitude = (atan2(unitVector.y(), unitVector.x()) * radians).wrap().asRadians();
103 } else {
104 _longitude = 0;
105 }
106}
107
108SpherePoint::SpherePoint() noexcept : _longitude(nan("")), _latitude(nan("")) {}
109
110SpherePoint::SpherePoint(SpherePoint const& other) noexcept = default;
111
112SpherePoint::SpherePoint(SpherePoint&& other) noexcept = default;
113
114SpherePoint& SpherePoint::operator=(SpherePoint const& other) noexcept = default;
115
116SpherePoint& SpherePoint::operator=(SpherePoint&& other) noexcept = default;
117
118SpherePoint::operator sphgeom::LonLat() const noexcept {
119 return sphgeom::LonLat::fromRadians(getLongitude().asRadians(), getLatitude().asRadians());
120}
121
122SpherePoint::~SpherePoint() noexcept = default;
123
124sphgeom::UnitVector3d SpherePoint::getVector() const noexcept {
125 return sphgeom::UnitVector3d::fromNormalized(cos(_longitude) * cos(_latitude),
126 sin(_longitude) * cos(_latitude), sin(_latitude));
127}
128
130 return Point2D(getLongitude().asAngularUnits(unit), getLatitude().asAngularUnits(unit));
131}
132
133Angle SpherePoint::operator[](size_t index) const {
134 switch (index) {
135 case 0:
136 return getLongitude();
137 case 1:
138 return getLatitude();
139 default:
140 throw pexExcept::OutOfRangeError("Index " + to_string(index) + " must be 0 or 1.");
141 }
142}
143
144bool SpherePoint::isFinite() const noexcept { return isfinite(_longitude) && isfinite(_latitude); }
145
146bool SpherePoint::operator==(SpherePoint const& other) const noexcept {
147 // Deliberate override of Style Guide 5-12
148 // Approximate FP comparison would make object equality intransitive
149 return _longitude == other._longitude && _latitude == other._latitude;
150}
151
152bool SpherePoint::operator!=(SpherePoint const& other) const noexcept { return !(*this == other); }
153
155 // Completely arbitrary seed
156 return utils::hashCombine(17, _longitude, _latitude);
157}
158
160 Angle const deltaLon = other.getLongitude() - this->getLongitude();
161
162 double const sinDelta1 = sin(getLatitude().asRadians());
163 double const sinDelta2 = sin(other.getLatitude().asRadians());
164 double const cosDelta1 = cos(getLatitude().asRadians());
165 double const cosDelta2 = cos(other.getLatitude().asRadians());
166
167 // Adapted from http://www.movable-type.co.uk/scripts/latlong.html
168 double const y = sin(deltaLon) * cosDelta2;
169 double const x = cosDelta1 * sinDelta2 - sinDelta1 * cosDelta2 * cos(deltaLon);
170 return (90.0 * degrees - atan2(y, x) * radians).wrap();
171}
172
173Angle SpherePoint::separation(SpherePoint const& other) const noexcept {
174 return haversine(getLongitude() - other.getLongitude(), getLatitude() - other.getLatitude(),
175 cos(getLatitude().asRadians()), cos(other.getLatitude().asRadians()));
176}
177
178SpherePoint SpherePoint::rotated(SpherePoint const& axis, Angle const& amount) const noexcept {
179 auto const rotation = Eigen::AngleAxisd(amount.asRadians(), asEigen(axis.getVector())).matrix();
180 auto const x = asEigen(getVector());
181 auto const xprime = rotation * x;
182 return SpherePoint(sphgeom::Vector3d(xprime[0], xprime[1], xprime[2]));
183}
184
185SpherePoint SpherePoint::offset(Angle const& bearing, Angle const& amount) const {
186 double const phi = bearing.asRadians();
187
188 // let v = vector in the direction bearing points (tangent to surface of sphere)
189 // To do the rotation, use rotated() method.
190 // - must provide an axis of rotation: take the cross product r x v to get that axis (pole)
191
192 auto r = getVector();
193
194 // Get the vector v:
195 // let u = unit vector lying on a parallel of declination
196 // let w = unit vector along line of longitude = r x u
197 // the vector v must satisfy the following:
198 // r . v = 0
199 // u . v = cos(phi)
200 // w . v = sin(phi)
201
202 // v is a linear combination of u and w
203 // v = cos(phi)*u + sin(phi)*w
204 sphgeom::Vector3d const u(-sin(_longitude), cos(_longitude), 0.0);
205 auto w = r.cross(u);
206 auto v = cos(phi) * u + sin(phi) * w;
207
208 // take r x v to get the axis
209 SpherePoint axis = SpherePoint(r.cross(v));
210
211 return rotated(axis, amount);
212}
213
215 geom::Angle const alpha1 = this->getLongitude();
216 geom::Angle const delta1 = this->getLatitude();
217 geom::Angle const alpha2 = other.getLongitude();
218 geom::Angle const delta2 = other.getLatitude();
219
220 // Compute the projection of "other" on a tangent plane centered at this point
221 double const sinDelta1 = std::sin(delta1);
222 double const cosDelta1 = std::cos(delta1);
223 double const sinDelta2 = std::sin(delta2);
224 double const cosDelta2 = std::cos(delta2);
225 double const cosAlphaDiff = std::cos(alpha2 - alpha1);
226 double const sinAlphaDiff = std::sin(alpha2 - alpha1);
227
228 double const div = cosDelta1 * cosAlphaDiff * cosDelta2 + sinDelta1 * sinDelta2;
229 double const xi = cosDelta1 * sinAlphaDiff / div;
230 double const eta = (cosDelta1 * cosAlphaDiff * sinDelta2 - sinDelta1 * cosDelta2) / div;
231
232 return std::make_pair(xi * geom::radians, eta * geom::radians);
233}
234
236 if (coords.size() == 0) {
237 throw LSST_EXCEPT(pex::exceptions::LengthError, "No coordinates provided to average");
238 }
239 sphgeom::Vector3d sum(0, 0, 0);
240 sphgeom::Vector3d corr(0, 0, 0); // Kahan summation correction
241 for (auto const& sp : coords) {
242 auto const point = sp.getVector();
243 // Kahan summation
244 auto const add = point - corr;
245 auto const temp = sum + add;
246 corr = (temp - sum) - add;
247 sum = temp;
248 }
249 sum /= static_cast<double>(coords.size());
250 return SpherePoint(sum);
251}
252
254 // Can't provide atomic guarantee anyway for I/O, so ok to set ostream states here.
255 auto oldFlags = os.setf(ostream::fixed);
256 // 10 digits of precision for a value in degrees is about 1 milliarcsecond
257 auto oldPrecision = os.precision(10);
258
259 os << "(" << point.getLongitude().asDegrees() << ", ";
260 os.setf(ostream::showpos);
261 os << point.getLatitude().asDegrees() << ")";
262
263 os.flags(oldFlags);
264 os.precision(oldPrecision);
265 return os;
266}
267
268} // namespace geom
269} // namespace lsst
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::ostream * os
Definition: Schema.cc:557
int y
Definition: SpanSet.cc:48
table::Key< lsst::geom::Angle > longitude
Definition: VisitInfo.cc:213
table::Key< lsst::geom::Angle > latitude
Definition: VisitInfo.cc:212
T asin(T... args)
T atan2(T... args)
A class representing an angle.
Definition: Angle.h:128
constexpr double asRadians() const noexcept
Return an Angle's value in radians.
Definition: Angle.h:173
constexpr double asDegrees() const noexcept
Return an Angle's value in degrees.
Definition: Angle.h:176
A class used to convert scalar POD types such as double to Angle.
Definition: Angle.h:71
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
SpherePoint & operator=(SpherePoint const &other) noexcept
Overwrite this object with the value of another SpherePoint.
std::pair< Angle, Angle > getTangentPlaneOffset(SpherePoint const &other) const
Get the offset from a tangent plane centered at this point to another point.
Definition: SpherePoint.cc:214
SpherePoint offset(Angle const &bearing, Angle const &amount) const
Return a point offset from this one along a great circle.
Definition: SpherePoint.cc:185
bool isFinite() const noexcept
true if this point is a well-defined position.
Definition: SpherePoint.cc:144
bool operator!=(SpherePoint const &other) const noexcept
false if two points represent the same position.
Definition: SpherePoint.cc:152
Point2D getPosition(AngleUnit unit) const noexcept
Return longitude, latitude as a Point2D object.
Definition: SpherePoint.cc:129
Angle bearingTo(SpherePoint const &other) const
Orientation at this point of the great circle arc to another point.
Definition: SpherePoint.cc:159
bool atPole() const noexcept
true if this point is either coordinate pole.
Definition: SpherePoint.h:237
Angle separation(SpherePoint const &other) const noexcept
Angular distance between two points.
Definition: SpherePoint.cc:173
sphgeom::UnitVector3d getVector() const noexcept
A unit vector representation of this point.
Definition: SpherePoint.cc:124
SpherePoint() noexcept
Construct a SpherePoint with "nan" for longitude and latitude.
Definition: SpherePoint.cc:108
Angle getLatitude() const noexcept
The latitude of this point.
Definition: SpherePoint.h:190
Angle getLongitude() const noexcept
The longitude of this point.
Definition: SpherePoint.h:178
bool operator==(SpherePoint const &other) const noexcept
true if two points have the same longitude and latitude.
Definition: SpherePoint.cc:146
SpherePoint rotated(SpherePoint const &axis, Angle const &amount) const noexcept
Return a point rotated from this one around an axis.
Definition: SpherePoint.cc:178
Angle operator[](size_t index) const
Longitude and latitude by index.
Definition: SpherePoint.cc:133
std::size_t hash_value() const noexcept
Return a hash of this object.
Definition: SpherePoint.cc:154
Reports invalid arguments.
Definition: Runtime.h:66
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
Reports attempts to access elements outside a valid range of indices.
Definition: Runtime.h:89
LonLat represents a spherical coordinate (longitude/latitude angle) pair.
Definition: LonLat.h:48
static LonLat fromRadians(double lon, double lat)
Definition: LonLat.h:55
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
static UnitVector3d fromNormalized(Vector3d const &v)
fromNormalized returns the unit vector equal to v, which is assumed to be normalized.
Definition: UnitVector3d.h:82
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition: Vector3d.h:44
T cos(T... args)
T div(T... args)
T fabs(T... args)
T flags(T... args)
T isfinite(T... args)
T make_pair(T... args)
T norm(const T &x)
Definition: Integrate.h:160
lsst::geom::Angle Angle
Definition: misc.h:33
lsst::geom::SpherePoint SpherePoint
Definition: misc.h:34
std::size_t hashCombine(std::size_t seed) noexcept
Combine hashes.
Definition: hashCombine.h:35
Eigen::Vector3d asEigen(sphgeom::Vector3d const &vector) noexcept
Definition: sphgeomUtils.h:36
SpherePoint averageSpherePoint(std::vector< SpherePoint > const &coords)
Return the average of a list of coordinates.
Definition: SpherePoint.cc:235
constexpr double HALFPI
Definition: Angle.h:42
constexpr AngleUnit degrees
constant with units of degrees
Definition: Angle.h:110
std::ostream & operator<<(std::ostream &os, lsst::geom::AffineTransform const &transform)
Point< double, 2 > Point2D
Definition: Point.h:324
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:109
def wrap(ctrl)
Definition: wrap.py:305
double sin(Angle const &a)
Definition: Angle.h:102
double cos(Angle const &a)
Definition: Angle.h:103
A base class for image defects.
STL namespace.
T precision(T... args)
T setf(T... args)
T sin(T... args)
T size(T... args)
T sqrt(T... args)
double w
Definition: CoaddPsf.cc:69
T str(T... args)