LSSTApplications  18.0.0+106,18.0.0+50,19.0.0,19.0.0+1,19.0.0+10,19.0.0+11,19.0.0+13,19.0.0+17,19.0.0+2,19.0.0-1-g20d9b18+6,19.0.0-1-g425ff20,19.0.0-1-g5549ca4,19.0.0-1-g580fafe+6,19.0.0-1-g6fe20d0+1,19.0.0-1-g7011481+9,19.0.0-1-g8c57eb9+6,19.0.0-1-gb5175dc+11,19.0.0-1-gdc0e4a7+9,19.0.0-1-ge272bc4+6,19.0.0-1-ge3aa853,19.0.0-10-g448f008b,19.0.0-12-g6990b2c,19.0.0-2-g0d9f9cd+11,19.0.0-2-g3d9e4fb2+11,19.0.0-2-g5037de4,19.0.0-2-gb96a1c4+3,19.0.0-2-gd955cfd+15,19.0.0-3-g2d13df8,19.0.0-3-g6f3c7dc,19.0.0-4-g725f80e+11,19.0.0-4-ga671dab3b+1,19.0.0-4-gad373c5+3,19.0.0-5-ga2acb9c+2,19.0.0-5-gfe96e6c+2,w.2020.01
LSSTDataManagementBasePackage
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 
28 #include "lsst/utils/hashCombine.h"
29 #include "lsst/geom/SpherePoint.h"
30 #include "lsst/geom/sphgeomUtils.h"
31 #include "lsst/pex/exceptions.h"
32 
34 using namespace std;
35 
36 namespace lsst {
37 
38 namespace geom {
40 namespace {
54 Angle 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.";
83  throw pexExcept::InvalidParameterError(buffer.str());
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 
97 void 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 
108 SpherePoint::SpherePoint() noexcept : _longitude(nan("")), _latitude(nan("")) {}
109 
110 SpherePoint::SpherePoint(SpherePoint const& other) noexcept = default;
111 
112 SpherePoint::SpherePoint(SpherePoint&& other) noexcept = default;
113 
114 SpherePoint& SpherePoint::operator=(SpherePoint const& other) noexcept = default;
115 
116 SpherePoint& SpherePoint::operator=(SpherePoint&& other) noexcept = default;
117 
118 SpherePoint::operator sphgeom::LonLat() const noexcept {
119  return sphgeom::LonLat::fromRadians(getLongitude().asRadians(), getLatitude().asRadians());
120 }
121 
122 SpherePoint::~SpherePoint() noexcept = default;
123 
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 
133 Angle 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 
144 bool SpherePoint::isFinite() const noexcept { return isfinite(_longitude) && isfinite(_latitude); }
145 
146 bool 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 
152 bool 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 
173 Angle 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 
178 SpherePoint 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 
185 SpherePoint 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 be sloppy.
255  auto oldFlags = os.setf(ostream::fixed);
256  auto oldPrecision = os.precision(6);
257 
258  os << "(" << point.getLongitude().asDegrees() << ", ";
259  os.setf(ostream::showpos);
260  os << point.getLatitude().asDegrees() << ")";
261 
262  os.flags(oldFlags);
263  os.precision(oldPrecision);
264  return os;
265 }
266 
267 } // namespace geom
268 } // namespace lsst
bool operator==(SpherePoint const &other) const noexcept
true if two points have the same longitude and latitude.
Definition: SpherePoint.cc:146
std::size_t hash_value() const noexcept
Return a hash of this object.
Definition: SpherePoint.cc:154
T setf(T... args)
bool isFinite() const noexcept
true if this point is a well-defined position.
Definition: SpherePoint.cc:144
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
table::Key< lsst::geom::Angle > latitude
Definition: VisitInfo.cc:171
Angle bearingTo(SpherePoint const &other) const
Orientation at this point of the great circle arc to another point.
Definition: SpherePoint.cc:159
constexpr double asRadians() const noexcept
Return an Angle&#39;s value in radians.
Definition: Angle.h:166
constexpr double asDegrees() const noexcept
Return an Angle&#39;s value in degrees.
Definition: Angle.h:169
Point2D getPosition(AngleUnit unit) const noexcept
Return longitude, latitude as a Point2D object.
Definition: SpherePoint.cc:129
Reports attempts to exceed implementation-defined length limits for some classes. ...
Definition: Runtime.h:76
int y
Definition: SpanSet.cc:49
T norm(const T &x)
Definition: Integrate.h:194
STL namespace.
T precision(T... args)
Angle getLongitude() const noexcept
The longitude of this point.
Definition: SpherePoint.h:178
SpherePoint rotated(SpherePoint const &axis, Angle const &amount) const noexcept
Return a point rotated from this one around an axis.
Definition: SpherePoint.cc:178
sphgeom::UnitVector3d getVector() const noexcept
A unit vector representation of this point.
Definition: SpherePoint.cc:124
ItemVariant const * other
Definition: Schema.cc:56
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition: Vector3d.h:44
double sin(Angle const &a)
Definition: Angle.h:102
AngleUnit constexpr radians
constant with units of radians
Definition: Angle.h:108
A class representing an angle.
Definition: Angle.h:127
double y() const
Definition: Vector3d.h:68
Angle operator[](size_t index) const
Longitude and latitude by index.
Definition: SpherePoint.cc:133
std::ostream & operator<<(std::ostream &os, SpherePoint const &point)
Print the value of a point to a stream.
Definition: SpherePoint.cc:253
Point< double, 2 > Point2D
Definition: Point.h:324
double cos(Angle const &a)
Definition: Angle.h:103
T atan2(T... args)
T sin(T... args)
T flags(T... args)
SpherePoint() noexcept
Construct a SpherePoint with "nan" for longitude and latitude.
Definition: SpherePoint.cc:108
double x() const
Definition: Vector3d.h:66
AngleUnit constexpr degrees
constant with units of degrees
Definition: Angle.h:109
A base class for image defects.
Angle separation(SpherePoint const &other) const noexcept
Angular distance between two points.
Definition: SpherePoint.cc:173
Angle getLatitude() const noexcept
The latitude of this point.
Definition: SpherePoint.h:190
T str(T... args)
T make_pair(T... args)
T isfinite(T... args)
SpherePoint offset(Angle const &bearing, Angle const &amount) const
Return a point offset from this one along a great circle.
Definition: SpherePoint.cc:185
T cos(T... args)
T fabs(T... args)
double x
T div(T... args)
Eigen::Vector3d asEigen(sphgeom::Vector3d const &vector) noexcept
Definition: sphgeomUtils.h:36
table::Key< lsst::geom::Angle > longitude
Definition: VisitInfo.cc:172
double w
Definition: CoaddPsf.cc:69
T asin(T... args)
T size(T... args)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Reports attempts to access elements outside a valid range of indices.
Definition: Runtime.h:89
STL class.
static LonLat fromRadians(double lon, double lat)
Definition: LonLat.h:55
SpherePoint & operator=(SpherePoint const &other) noexcept
Overwrite this object with the value of another SpherePoint.
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
LonLat represents a spherical coordinate (longitude/latitude angle) pair.
Definition: LonLat.h:48
double constexpr HALFPI
Definition: Angle.h:41
bool atPole() const noexcept
true if this point is either coordinate pole.
Definition: SpherePoint.h:237
def wrap(ctrl)
Definition: wrap.py:302
double z() const
Definition: Vector3d.h:70
Reports invalid arguments.
Definition: Runtime.h:66
lsst::geom::SpherePoint SpherePoint
Definition: misc.h:35
lsst::geom::Angle Angle
Definition: misc.h:33
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
bool operator!=(SpherePoint const &other) const noexcept
false if two points represent the same position.
Definition: SpherePoint.cc:152
T sqrt(T... args)
std::size_t hashCombine(std::size_t seed) noexcept
Combine hashes.
Definition: hashCombine.h:35
std::ostream * os
Definition: Schema.cc:746
STL class.
static UnitVector3d fromNormalized(Vector3d const &v)
fromNormalized returns the unit vector equal to v, which is assumed to be normalized.
Definition: UnitVector3d.h:82
A class used to convert scalar POD types such as double to Angle.
Definition: Angle.h:70
SpherePoint averageSpherePoint(std::vector< SpherePoint > const &coords)
Return the average of a list of coordinates.
Definition: SpherePoint.cc:235
double getNorm() const
getNorm returns the L2 norm of this vector.
Definition: Vector3d.h:81