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