LSST Applications g1cfbe01172+01aa18f939,g20cdd03214+31e6b93548,g28da252d5a+ea8665a95b,g2bbee38e9b+9ec6cc348d,g2bc492864f+9ec6cc348d,g347aa1857d+9ec6cc348d,g3a166c0a6a+9ec6cc348d,g4322eb9e3a+65eff1e020,g461a3dce89+b86e4b8053,g50ff169b8f+f991eae79d,g52b1c1532d+b86e4b8053,g607f77f49a+31e6b93548,g78056777b3+8ae2798781,g858d7b2824+31e6b93548,g8cd86fa7b1+4851e61ca4,g9ddcbc5298+f24b38b85a,ga1e77700b3+3309dba821,gae0086650b+b86e4b8053,gb0e22166c9+6076c0b52b,gbb886bcc26+dccb771098,gbd462c55f0+dc07f8e65d,gc0c51c7ec2+31e6b93548,gc120e1dc64+a417ce3171,gc28159a63d+9ec6cc348d,gc2a6998b7e+f95f64aeae,gcdd4ae20e8+507450c4cd,gcf0d15dbbd+507450c4cd,gd1535ee943+bcf88ba65f,gd598c5cd71+66126f91fb,gdaeeff99f8+006e14e809,gdbce86181e+39d5515b1a,ge3d4d395c2+b12d4d6a95,ge79ae78c31+9ec6cc348d,gf048a9a2f4+d9c36e6b63,gfbcc870c63+ea41c4420b,w.2024.27
LSST Data Management Base Package
Loading...
Searching...
No Matches
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 =
90 sphgeom::UnitVector3d::fromNormalized(vector.x() / norm, vector.y() / norm, vector.z() / norm);
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
#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:215
table::Key< lsst::geom::Angle > latitude
Definition VisitInfo.cc:214
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.
SpherePoint offset(Angle const &bearing, Angle const &amount) const
Return a point offset from this one along a great circle.
bool isFinite() const noexcept
true if this point is a well-defined position.
bool operator!=(SpherePoint const &other) const noexcept
false if two points represent the same position.
Point2D getPosition(AngleUnit unit) const noexcept
Return longitude, latitude as a Point2D object.
Angle bearingTo(SpherePoint const &other) const
Orientation at this point of the great circle arc to another point.
bool atPole() const noexcept
true if this point is either coordinate pole.
Angle separation(SpherePoint const &other) const noexcept
Angular distance between two points.
sphgeom::UnitVector3d getVector() const noexcept
A unit vector representation of this point.
SpherePoint() noexcept
Construct a SpherePoint with "nan" for longitude and latitude.
Angle getLatitude() const noexcept
The latitude of this point.
Angle getLongitude() const noexcept
The longitude of this point.
bool operator==(SpherePoint const &other) const noexcept
true if two points have the same longitude and latitude.
SpherePoint rotated(SpherePoint const &axis, Angle const &amount) const noexcept
Return a point rotated from this one around an axis.
Angle operator[](size_t index) const
Longitude and latitude by index.
std::size_t hash_value() const noexcept
Return a hash of this object.
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:55
static LonLat fromRadians(double lon, double lat)
Definition LonLat.h:62
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
static UnitVector3d fromNormalized(Vector3d const &v)
fromNormalized returns the unit vector equal to v, which is assumed to be normalized.
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition Vector3d.h:51
T cos(T... args)
T div(T... args)
T fabs(T... args)
T isfinite(T... args)
T make_pair(T... args)
std::size_t hashCombine(std::size_t seed) noexcept
Combine hashes.
Definition hashCombine.h:35
Eigen::Vector3d asEigen(sphgeom::Vector3d const &vector) noexcept
double constexpr HALFPI
Definition Angle.h:42
SpherePoint averageSpherePoint(std::vector< SpherePoint > const &coords)
Return the average of a list of coordinates.
AngleUnit constexpr radians
constant with units of radians
Definition Angle.h:109
std::ostream & operator<<(std::ostream &os, lsst::geom::AffineTransform const &transform)
Point< double, 2 > Point2D
Definition Point.h:324
STL namespace.
T nan(T... args)
T setf(T... args)
T sin(T... args)
T size(T... args)
T sqrt(T... args)
double w
Definition CoaddPsf.cc:70
T str(T... args)
T to_string(T... args)