LSST Applications g063fba187b+66a50001ff,g0f08755f38+1a22dc2551,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g217e2c1bcf+12e87a5bd8,g246886dfd9+466c7b6c06,g28da252d5a+858b171e04,g2bbee38e9b+404b60ec9b,g2bc492864f+404b60ec9b,g3156d2b45e+6e55a43351,g347aa1857d+404b60ec9b,g35bb328faa+a8ce1bb630,g3a166c0a6a+404b60ec9b,g3e281a1b8c+c5dd892a6c,g414038480c+6b9177ef31,g41af890bb2+70bea58702,g599934f4f4+b8c5400ca5,g781aacb6e4+a8ce1bb630,g7af13505b9+b5b9cefdb8,g80478fca09+c2997882f3,g82479be7b0+8974e6af0f,g858d7b2824+1a22dc2551,g89c8672015+f4add4ffd5,g8f1c07a47a+de51c9b0a5,g9125e01d80+a8ce1bb630,ga5288a1d22+b66f8cf76b,gb58c049af0+d64f4d3760,gc28159a63d+404b60ec9b,gcab2d0539d+66cf1de5d4,gcf0d15dbbd+12cb7e2563,gda6a2b7d83+12cb7e2563,gdaeeff99f8+1711a396fd,ge79ae78c31+404b60ec9b,gef2f8181fd+414189b318,gf0baf85859+c1f95f4921,gf0c06eb49c+1a22dc2551,gfa517265be+1a22dc2551,gfa999e8aa5+17cd334064,v28.0.0.rc2
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 cpputils::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)