LSSTApplications  19.0.0-14-gb0260a2+72efe9b372,20.0.0+7927753e06,20.0.0+8829bf0056,20.0.0+995114c5d2,20.0.0+b6f4b2abd1,20.0.0+bddc4f4cbe,20.0.0-1-g253301a+8829bf0056,20.0.0-1-g2b7511a+0d71a2d77f,20.0.0-1-g5b95a8c+7461dd0434,20.0.0-12-g321c96ea+23efe4bbff,20.0.0-16-gfab17e72e+fdf35455f6,20.0.0-2-g0070d88+ba3ffc8f0b,20.0.0-2-g4dae9ad+ee58a624b3,20.0.0-2-g61b8584+5d3db074ba,20.0.0-2-gb780d76+d529cf1a41,20.0.0-2-ged6426c+226a441f5f,20.0.0-2-gf072044+8829bf0056,20.0.0-2-gf1f7952+ee58a624b3,20.0.0-20-geae50cf+e37fec0aee,20.0.0-25-g3dcad98+544a109665,20.0.0-25-g5eafb0f+ee58a624b3,20.0.0-27-g64178ef+f1f297b00a,20.0.0-3-g4cc78c6+e0676b0dc8,20.0.0-3-g8f21e14+4fd2c12c9a,20.0.0-3-gbd60e8c+187b78b4b8,20.0.0-3-gbecbe05+48431fa087,20.0.0-38-ge4adf513+a12e1f8e37,20.0.0-4-g97dc21a+544a109665,20.0.0-4-gb4befbc+087873070b,20.0.0-4-gf910f65+5d3db074ba,20.0.0-5-gdfe0fee+199202a608,20.0.0-5-gfbfe500+d529cf1a41,20.0.0-6-g64f541c+d529cf1a41,20.0.0-6-g9a5b7a1+a1cd37312e,20.0.0-68-ga3f3dda+5fca18c6a4,20.0.0-9-g4aef684+e18322736b,w.2020.45
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 =
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 
124 sphgeom::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 
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 
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
y
int y
Definition: SpanSet.cc:49
lsst::geom::degrees
constexpr AngleUnit degrees
constant with units of degrees
Definition: Angle.h:109
lsst::geom::averageSpherePoint
SpherePoint averageSpherePoint(std::vector< SpherePoint > const &coords)
Return the average of a list of coordinates.
Definition: SpherePoint.cc:235
lsst::geom::SpherePoint::getPosition
Point2D getPosition(AngleUnit unit) const noexcept
Return longitude, latitude as a Point2D object.
Definition: SpherePoint.cc:129
lsst::sphgeom::sin
double sin(Angle const &a)
Definition: Angle.h:102
sphgeomUtils.h
std::fabs
T fabs(T... args)
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:194
std::atan2
T atan2(T... args)
lsst::afw::table::SpherePoint
lsst::geom::SpherePoint SpherePoint
Definition: misc.h:35
std::pair
std::cos
T cos(T... args)
std::asin
T asin(T... args)
lsst::sphgeom::UnitVector3d::x
double x() const
Definition: UnitVector3d.h:145
lsst::geom::SpherePoint::hash_value
std::size_t hash_value() const noexcept
Return a hash of this object.
Definition: SpherePoint.cc:154
std::vector
STL class.
std::vector::size
T size(T... args)
lsst::sphgeom::LonLat
LonLat represents a spherical coordinate (longitude/latitude angle) pair.
Definition: LonLat.h:48
lsst::geom::SpherePoint::separation
Angle separation(SpherePoint const &other) const noexcept
Angular distance between two points.
Definition: SpherePoint.cc:173
lsst::geom::SpherePoint::SpherePoint
SpherePoint() noexcept
Construct a SpherePoint with "nan" for longitude and latitude.
Definition: SpherePoint.cc:108
std::div
T div(T... args)
lsst::geom::SpherePoint::getLongitude
Angle getLongitude() const noexcept
The longitude of this point.
Definition: SpherePoint.h:178
std::stringstream
STL class.
lsst::geom::SpherePoint::operator==
bool operator==(SpherePoint const &other) const noexcept
true if two points have the same longitude and latitude.
Definition: SpherePoint.cc:146
lsst::geom::SpherePoint::getTangentPlaneOffset
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
lsst::geom::Point2D
Point< double, 2 > Point2D
Definition: Point.h:324
std::sqrt
T sqrt(T... args)
lsst::sphgeom::Vector3d
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition: Vector3d.h:44
lsst::geom::SpherePoint::operator[]
Angle operator[](size_t index) const
Longitude and latitude by index.
Definition: SpherePoint.cc:133
lsst::geom::Angle::asRadians
constexpr double asRadians() const noexcept
Return an Angle's value in radians.
Definition: Angle.h:166
hashCombine.h
lsst::afw::table::Angle
lsst::geom::Angle Angle
Definition: misc.h:33
lsst::sphgeom::UnitVector3d
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
lsst::geom::SpherePoint::getVector
sphgeom::UnitVector3d getVector() const noexcept
A unit vector representation of this point.
Definition: SpherePoint.cc:124
lsst::geom::SpherePoint::bearingTo
Angle bearingTo(SpherePoint const &other) const
Orientation at this point of the great circle arc to another point.
Definition: SpherePoint.cc:159
std::isfinite
T isfinite(T... args)
std::ostream
STL class.
x
double x
Definition: ChebyshevBoundedField.cc:277
lsst::geom::SpherePoint::rotated
SpherePoint rotated(SpherePoint const &axis, Angle const &amount) const noexcept
Return a point rotated from this one around an axis.
Definition: SpherePoint.cc:178
lsst.pex::exceptions::LengthError
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
other
ItemVariant const * other
Definition: Schema.cc:56
std::ostream::setf
T setf(T... args)
lsst::geom::SpherePoint::operator=
SpherePoint & operator=(SpherePoint const &other) noexcept
Overwrite this object with the value of another SpherePoint.
std::ostream::flags
T flags(T... args)
lsst::geom::SpherePoint::atPole
bool atPole() const noexcept
true if this point is either coordinate pole.
Definition: SpherePoint.h:237
lsst::sphgeom::UnitVector3d::z
double z() const
Definition: UnitVector3d.h:149
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
lsst::geom::operator<<
std::ostream & operator<<(std::ostream &os, lsst::geom::AffineTransform const &transform)
Definition: AffineTransform.cc:72
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::sin
T sin(T... args)
lsst::geom
Definition: AffineTransform.h:36
os
std::ostream * os
Definition: Schema.cc:746
lsst::utils::hashCombine
std::size_t hashCombine(std::size_t seed) noexcept
Combine hashes.
Definition: hashCombine.h:35
lsst::geom::SpherePoint::operator!=
bool operator!=(SpherePoint const &other) const noexcept
false if two points represent the same position.
Definition: SpherePoint.cc:152
lsst.pex::exceptions::InvalidParameterError
Reports invalid arguments.
Definition: Runtime.h:66
lsst::sphgeom::UnitVector3d::y
double y() const
Definition: UnitVector3d.h:147
std
STL namespace.
lsst::geom::Point< double, 2 >
lsst.pex.config.wrap.wrap
def wrap(ctrl)
Definition: wrap.py:302
lsst::geom::AngleUnit
A class used to convert scalar POD types such as double to Angle.
Definition: Angle.h:70
lsst.pex::exceptions::OutOfRangeError
Reports attempts to access elements outside a valid range of indices.
Definition: Runtime.h:89
lsst::geom::Angle
A class representing an angle.
Definition: Angle.h:127
lsst::geom::SpherePoint::getLatitude
Angle getLatitude() const noexcept
The latitude of this point.
Definition: SpherePoint.h:190
lsst::geom::radians
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:108
lsst::geom::SpherePoint::isFinite
bool isFinite() const noexcept
true if this point is a well-defined position.
Definition: SpherePoint.cc:144
lsst::geom::asEigen
Eigen::Vector3d asEigen(sphgeom::Vector3d const &vector) noexcept
Definition: sphgeomUtils.h:36
lsst::geom::Angle::asDegrees
constexpr double asDegrees() const noexcept
Return an Angle's value in degrees.
Definition: Angle.h:169
std::stringstream::str
T str(T... args)
w
double w
Definition: CoaddPsf.cc:69
lsst.pex::exceptions
Definition: Exception.h:37
std::size_t
SpherePoint.h
std::make_pair
T make_pair(T... args)
longitude
table::Key< lsst::geom::Angle > longitude
Definition: VisitInfo.cc:172
lsst::geom::SpherePoint
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
latitude
table::Key< lsst::geom::Angle > latitude
Definition: VisitInfo.cc:171
astshim.fitsChanContinued.to_string
def to_string(self)
Definition: fitsChanContinued.py:117
lsst::sphgeom::LonLat::fromRadians
static LonLat fromRadians(double lon, double lat)
Definition: LonLat.h:55
lsst::geom::SpherePoint::offset
SpherePoint offset(Angle const &bearing, Angle const &amount) const
Return a point offset from this one along a great circle.
Definition: SpherePoint.cc:185
lsst::geom::SpherePoint::~SpherePoint
~SpherePoint() noexcept
lsst::sphgeom::UnitVector3d::fromNormalized
static UnitVector3d fromNormalized(Vector3d const &v)
fromNormalized returns the unit vector equal to v, which is assumed to be normalized.
Definition: UnitVector3d.h:82
lsst::sphgeom::cos
double cos(Angle const &a)
Definition: Angle.h:103
lsst::geom::HALFPI
constexpr double HALFPI
Definition: Angle.h:41
std::ostream::precision
T precision(T... args)
exceptions.h