SLongitudeLatitude.h
Go to the documentation of this file.
1 
2 #pragma once
3 
4 #include <ionMath.h>
5 
6 
7 namespace ion
8 {
9 
10  enum class ECompassDirection
11  {
12  N = 1, E = 1,
13  S = -1, W = -1
14  };
15 
16  template <typename T>
17  struct SLongitudeLatitude : public vec2<T>
18  {
19 
22 
25  : Longitude(this->X), Latitude(this->Y), vec2<T>()
26  {}
27 
29  SLongitudeLatitude(T const lon, T const lat)
30  : Longitude(this->X), Latitude(this->Y), vec2<T>(lon, lat)
31  {}
32 
35  : Longitude(this->X), Latitude(this->Y), vec2<T>(vec)
36  {}
37 
39  SLongitudeLatitude<T> & operator = (SLongitudeLatitude<T> const & vec)
40  {
41  this->Longitude = vec.Longitude;
42  this->Latitude = vec.Latitude;
43 
44  return * this;
45  }
46 
47  static T DMStoDecimal(T const Degrees, T const Minutes, T const Seconds, ECompassDirection const Direction = ECompassDirection::N)
48  {
49  return (Degrees + Minutes / 60 + Seconds / 3600) * (int) Direction;
50  }
51 
52  static void DecimalToDMS(T const Decimal, T & Degrees, T & Minutes, T & Seconds)
53  {
54  T Remainder = Decimal;
55 
56  Degrees = floor(Remainder);
57  Remainder -= Degrees;
58 
59  Remainder *= 60;
60  Minutes = floor(Remainder);
61  Remainder -= Minutes;
62 
63  Remainder *= 60;
64  Seconds = Remainder;
65  }
66 
67  static T DMStoDecimal(std::string const & String)
68  {
69  double Deg, Min, Sec;
70  char Dir;
71  ECompassDirection Direction;
72 
73  sscanf(String.c_str(), "%lf %lf %lf %c", & Deg, & Min, & Sec, & Dir);
74 
75  switch (tolower(Dir))
76  {
77  default:
78  case 'N':
79  Direction = ECompassDirection::N;
80  break;
81  case 'E':
82  Direction = ECompassDirection::E;
83  break;
84  case 'S':
85  Direction = ECompassDirection::S;
86  break;
87  case 'W':
88  Direction = ECompassDirection::W;
89  break;
90  }
91 
92  return DMStoDecimal(Deg, Min, Sec, Direction);
93  }
94 
95  vec2f Vector() const
96  {
97  return vec2f(Longitude, Latitude);
98  }
99 
101  {
102 
103  public:
104 
105  virtual T DistanceBetween(SLongitudeLatitude const &, SLongitudeLatitude const &) = 0;
106  virtual vec2<T> OffsetBetween(SLongitudeLatitude const &, SLongitudeLatitude const &) = 0;
107 
108  };
109 
111  {
112 
113  public:
114 
115  enum class EOffsetMode
116  {
117  Left,
118  Right,
119  Average
120  };
121 
123  {
124  T EarthRadius = (T) 6378.137;
125  T DeltaLat = DegreesToRadians(Right.Latitude - Left.Latitude);
126  T DeltaLong = DegreesToRadians(Right.Longitude - Left.Longitude);
127  T A =
128  Sin(DeltaLat / 2) * Sin(DeltaLat / 2) +
129  Cos(DegreesToRadians(Left.Latitude)) * Cos(DegreesToRadians(Right.Latitude)) * Sin(DeltaLong / 2) * Sin(DeltaLong / 2);
130  T C = 2 * ArcTan(Sqrt(A), Sqrt(1 - A));
131  T Distance = EarthRadius * C;
132 
133  return Distance * 1000;
134  }
135 
137  {
138  vec2<T> Offset;
140 
141  Left.Longitude = A.Longitude;
142  Right.Longitude = B.Longitude;
143  switch (OffsetMode)
144  {
145  default:
147  Left.Latitude = Right.Latitude = (A.Latitude + B.Latitude) / 2;
148  break;
149  case EOffsetMode::Left:
150  Left.Latitude = Right.Latitude = A.Latitude;
151  break;
152  case EOffsetMode::Right:
153  Left.Latitude = Right.Latitude = B.Latitude;
154  break;
155  }
156  Offset.X = (Left.Longitude < Right.Longitude ? 1 : -1) * DistanceBetween(Left, Right);
157 
158  Left.Latitude = A.Latitude;
159  Right.Latitude = B.Latitude;
160  switch (OffsetMode)
161  {
162  default:
164  Left.Longitude = Right.Longitude = (A.Longitude + B.Longitude) / 2;
165  break;
166  case EOffsetMode::Left:
167  Left.Longitude = Right.Longitude = A.Longitude;
168  break;
169  case EOffsetMode::Right:
170  Left.Longitude = Right.Longitude = B.Longitude;
171  break;
172  }
173  Offset.Y = (Left.Latitude < Right.Latitude ? 1 : -1) * DistanceBetween(Left, Right);
174 
175  return Offset;
176  }
177 
179  : OffsetMode(EOffsetMode::Average)
180  {}
181 
183  : OffsetMode(OffsetMode)
184  {}
185 
187 
188  };
189 
191  {
192 
193  public:
194 
196  {
197  T const lon1 = DegreesToRadians(Left.Longitude), lat1 = DegreesToRadians(Left.Latitude);
198  T const lon2 = DegreesToRadians(Right.Longitude), lat2 = DegreesToRadians(Right.Latitude);
199 
200  T const a = (T) 6378137, b = (T) 6356752.314245, f = (T) (1 / 298.257223563); // WGS-84 ellipsoid params
201  T const L = lon2 - lon1;
202  T const U1 = ArcTan((1 - f) * Tan(lat1));
203  T const U2 = ArcTan((1 - f) * Tan(lat2));
204  T const sinU1 = Sin(U1), cosU1 = Cos(U1);
205  T const sinU2 = Sin(U2), cosU2 = Cos(U2);
206 
207  T lambda = L, lambdaP, iterLimit = 100;
208  T cosSqAlpha, sinSigma, cos2SigmaM, cosSigma, sigma;
209  do {
210  T sinLambda = Sin(lambda), cosLambda = Cos(lambda);
211  sinSigma = Sqrt((cosU2*sinLambda) * (cosU2*sinLambda) +
212  (cosU1*sinU2 - sinU1*cosU2*cosLambda) * (cosU1*sinU2 - sinU1*cosU2*cosLambda));
213  if (sinSigma == 0)
214  return 0; // co-incident points
215  cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
216  sigma = ArcTan(sinSigma, cosSigma);
217  T sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
218  cosSqAlpha = 1 - sinAlpha*sinAlpha;
219  cos2SigmaM = cosSigma - 2 * sinU1*sinU2 / cosSqAlpha;
220  if (IsNaN(cos2SigmaM))
221  cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (6)
222  T C = f / 16 * cosSqAlpha*(4 + f*(4 - 3 * cosSqAlpha));
223  lambdaP = lambda;
224  lambda = L + (1 - C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM + C*cosSigma*(-1 + 2 * cos2SigmaM*cos2SigmaM)));
225  } while (Abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0);
226 
227  if (iterLimit == 0)
228  {
229  std::cerr << "Failed to converge in Vincenty Long/Lat conversion." << std::endl;
230  return std::numeric_limits<T>::max(); // formula failed to converge
231  }
232 
233  T uSq = cosSqAlpha * (a*a - b*b) / (b*b);
234  T A = 1 + uSq / 16384 * (4096 + uSq*(-768 + uSq*(320 - 175 * uSq)));
235  T B = uSq / 1024 * (256 + uSq*(-128 + uSq*(74 - 47 * uSq)));
236  T deltaSigma = B*sinSigma*(cos2SigmaM + B / 4 * (cosSigma*(-1 + 2 * cos2SigmaM*cos2SigmaM) - B / 6 * cos2SigmaM*(-3 + 4 * sinSigma*sinSigma)*(-3 + 4 * cos2SigmaM*cos2SigmaM)));
237  T s = b*A*(sigma - deltaSigma);
238 
239  return s;
240  }
241 
243  {}
244 
245  // Still causing problems on linux...
246  //CVincentyProjection(EOffsetMode const offsetMode)
247  // : CHaversineProjection(OffsetMode)
248  //{}
249 
250  };
251 
253  {
254 
255  public:
256 
258  {
259  CHaversineProjection Haversine;
260  SLongitudeLatitude const Center((Left.Longitude + Right.Longitude) / 2, PhiStandard);
261  vec2<T> const UnitLength = Haversine.OffsetBetween(Center - 0.5, Center + 0.5);
262 
263  vec2<T> LeftProjected(Left.Longitude * Cos(PhiStandard), Left.Latitude);
264  vec2<T> RightProjected(Right.Longitude * Cos(PhiStandard), Right.Latitude);
265 
266  LeftProjected *= UnitLength;
267  RightProjected *= UnitLength;
268 
269  return LeftProjected.GetDistanceFrom(RightProjected);
270  }
271 
273  {
274  CHaversineProjection Haversine;
275  SLongitudeLatitude const Center((Left.Longitude + Right.Longitude) / 2, PhiStandard);
276  vec2<T> const UnitLength = Haversine.OffsetBetween(Center - 0.5, Center + 0.5);
277 
278  vec2<T> LeftProjected(Left.Longitude * Cos(PhiStandard), Left.Latitude);
279  vec2<T> RightProjected(Right.Longitude * Cos(PhiStandard), Right.Latitude);
280 
281  LeftProjected *= UnitLength;
282  RightProjected *= UnitLength;
283 
284  return RightProjected - LeftProjected;
285  }
286 
287  CEquirectangularProjection(T const phiStandard)
288  : PhiStandard(phiStandard)
289  {}
290 
292 
293  };
294 
295  T DistanceTo(SLongitudeLatitude const & Other, SharedPointer<IProjectionSystem> Projection = SharedFromNew(new CHaversineProjection())) const
296  {
297  return Projection->DistanceBetween(* this, Other);
298  }
299 
300  vec2<T> OffsetTo(SLongitudeLatitude const & Other, SharedPointer<IProjectionSystem> Projection = SharedFromNew(new CHaversineProjection())) const
301  {
302  return Projection->OffsetBetween(* this, Other);
303  }
304 
305  };
306 
307 
310 
311 }
T GetDistanceFrom(vec2< T > const &v) const
Definition: vec2.h:179
T Abs(T const value)
Definition: ionMathFunctions.h:13
2D vector
Definition: vec2.h:17
static T DMStoDecimal(T const Degrees, T const Minutes, T const Seconds, ECompassDirection const Direction=ECompassDirection::N)
Definition: SLongitudeLatitude.h:47
T Average(T const a, T const b)
Definition: ionUtils.h:59
static void DecimalToDMS(T const Decimal, T &Degrees, T &Minutes, T &Seconds)
Definition: SLongitudeLatitude.h:52
float const e
Definition: ionConstants.h:15
SLongitudeLatitude< double > longlatd
Definition: SLongitudeLatitude.h:309
static T DMStoDecimal(std::string const &String)
Definition: SLongitudeLatitude.h:67
SLongitudeLatitude()
Default constructor.
Definition: SLongitudeLatitude.h:24
virtual T DistanceBetween(SLongitudeLatitude const &Left, SLongitudeLatitude const &Right)
Definition: SLongitudeLatitude.h:257
T DistanceTo(SLongitudeLatitude const &Other, SharedPointer< IProjectionSystem > Projection=SharedFromNew(new CHaversineProjection())) const
Definition: SLongitudeLatitude.h:295
vec2< T > OffsetTo(SLongitudeLatitude const &Other, SharedPointer< IProjectionSystem > Projection=SharedFromNew(new CHaversineProjection())) const
Definition: SLongitudeLatitude.h:300
T X
Definition: vec2.h:20
vec2< float > vec2f
Definition: vec2.h:574
Definition: CCatmullRomAdvancedSplineInterpolator.h:7
T & Latitude
Definition: SLongitudeLatitude.h:21
ECompassDirection
Definition: SLongitudeLatitude.h:10
virtual T DistanceBetween(SLongitudeLatitude const &Left, SLongitudeLatitude const &Right)
Definition: SLongitudeLatitude.h:122
T ArcTan(T const value)
Definition: ionMathFunctions.h:66
virtual vec2< T > OffsetBetween(SLongitudeLatitude const &A, SLongitudeLatitude const &B)
Definition: SLongitudeLatitude.h:136
SLongitudeLatitude< float > longlatf
Definition: SLongitudeLatitude.h:308
CHaversineProjection()
Definition: SLongitudeLatitude.h:178
CHaversineProjection(EOffsetMode const offsetMode)
Definition: SLongitudeLatitude.h:182
T & Longitude
Definition: SLongitudeLatitude.h:20
Definition: SLongitudeLatitude.h:110
Helper methods for dealing with strings.
Definition: String.h:15
EOffsetMode
Definition: SLongitudeLatitude.h:115
Definition: SLongitudeLatitude.h:100
SLongitudeLatitude(SLongitudeLatitude< T > const &vec)
Copy constructor.
Definition: SLongitudeLatitude.h:34
Definition: SLongitudeLatitude.h:252
CVincentyProjection()
Definition: SLongitudeLatitude.h:242
virtual vec2< T > OffsetBetween(SLongitudeLatitude const &Left, SLongitudeLatitude const &Right)
Definition: SLongitudeLatitude.h:272
T Tan(T const value)
Definition: ionMathFunctions.h:59
T Y
Definition: vec2.h:21
bool IsNaN(T const value)
Definition: ionMathFunctions.h:101
T Sin(T const value)
Definition: ionMathFunctions.h:31
vec2f Vector() const
Definition: SLongitudeLatitude.h:95
Definition: SLongitudeLatitude.h:190
SLongitudeLatitude(T const lon, T const lat)
Explicit constructor.
Definition: SLongitudeLatitude.h:29
EOffsetMode OffsetMode
Definition: SLongitudeLatitude.h:186
T Cos(T const value)
Definition: ionMathFunctions.h:45
Definition: SLongitudeLatitude.h:17
std::shared_ptr< T > SharedPointer
Definition: ionSmartPointer.h:25
T PhiStandard
Definition: SLongitudeLatitude.h:291
T Sqrt(T const value)
Definition: ionMathFunctions.h:80
CEquirectangularProjection(T const phiStandard)
Definition: SLongitudeLatitude.h:287
virtual T DistanceBetween(SLongitudeLatitude const &Left, SLongitudeLatitude const &Right)
Definition: SLongitudeLatitude.h:195