Subversion Repositories Projects

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2287 - 1

2
namespace GMap.NET
3
{
4
   using System;
5
   using System.Collections.Generic;
6
   using System.Diagnostics;
7
 
8
   /// <summary>
9
   /// defines projection
10
   /// </summary>
11
   public abstract class PureProjection
12
   {
13
      readonly List<Dictionary<PointLatLng, GPoint>> FromLatLngToPixelCache = new List<Dictionary<PointLatLng, GPoint>>(33);
14
      readonly List<Dictionary<GPoint, PointLatLng>> FromPixelToLatLngCache = new List<Dictionary<GPoint, PointLatLng>>(33);
15
 
16
      public PureProjection()
17
      {
18
         for(int i = 0; i < FromLatLngToPixelCache.Capacity; i++)
19
         {
20
            FromLatLngToPixelCache.Add(new Dictionary<PointLatLng, GPoint>());
21
            FromPixelToLatLngCache.Add(new Dictionary<GPoint, PointLatLng>());
22
         }
23
      }
24
 
25
      /// <summary>
26
      /// size of tile
27
      /// </summary>
28
      public abstract GSize TileSize
29
      {
30
         get;
31
      }
32
 
33
      /// <summary>
34
      /// Semi-major axis of ellipsoid, in meters
35
      /// </summary>
36
      public abstract double Axis
37
      {
38
         get;
39
      }
40
 
41
      /// <summary>
42
      /// Flattening of ellipsoid
43
      /// </summary>
44
      public abstract double Flattening
45
      {
46
         get;
47
      }
48
 
49
      /// <summary>
50
      /// get pixel coordinates from lat/lng
51
      /// </summary>
52
      /// <param name="lat"></param>
53
      /// <param name="lng"></param>
54
      /// <param name="zoom"></param>
55
      /// <returns></returns>
56
      public abstract GPoint FromLatLngToPixel(double lat, double lng, int zoom);
57
 
58
      /// <summary>
59
      /// gets lat/lng coordinates from pixel coordinates
60
      /// </summary>
61
      /// <param name="x"></param>
62
      /// <param name="y"></param>
63
      /// <param name="zoom"></param>
64
      /// <returns></returns>
65
      public abstract PointLatLng FromPixelToLatLng(long x, long y, int zoom);
66
 
67
      public GPoint FromLatLngToPixel(PointLatLng p, int zoom)
68
      {
69
         return FromLatLngToPixel(p, zoom, false);
70
      }
71
 
72
      /// <summary>
73
      /// get pixel coordinates from lat/lng
74
      /// </summary>
75
      /// <param name="p"></param>
76
      /// <param name="zoom"></param>
77
      /// <returns></returns>
78
      public GPoint FromLatLngToPixel(PointLatLng p, int zoom, bool useCache)
79
      {
80
         if(useCache)
81
         {
82
            GPoint ret = GPoint.Empty;
83
            if(!FromLatLngToPixelCache[zoom].TryGetValue(p, out ret))
84
            {
85
               ret = FromLatLngToPixel(p.Lat, p.Lng, zoom);
86
               FromLatLngToPixelCache[zoom].Add(p, ret);
87
 
88
               // for reverse cache
89
               if(!FromPixelToLatLngCache[zoom].ContainsKey(ret))
90
               {
91
                  FromPixelToLatLngCache[zoom].Add(ret, p);
92
               }
93
 
94
               Debug.WriteLine("FromLatLngToPixelCache[" + zoom + "] added " + p + " with " + ret);
95
            }
96
            return ret;
97
         }
98
         else
99
         {
100
            return FromLatLngToPixel(p.Lat, p.Lng, zoom);
101
         }
102
      }
103
 
104
      public PointLatLng FromPixelToLatLng(GPoint p, int zoom)
105
      {
106
         return FromPixelToLatLng(p, zoom, false);
107
      }
108
 
109
      /// <summary>
110
      /// gets lat/lng coordinates from pixel coordinates
111
      /// </summary>
112
      /// <param name="p"></param>
113
      /// <param name="zoom"></param>
114
      /// <returns></returns>
115
      public PointLatLng FromPixelToLatLng(GPoint p, int zoom, bool useCache)
116
      {
117
         if(useCache)
118
         {
119
            PointLatLng ret = PointLatLng.Empty;
120
            if(!FromPixelToLatLngCache[zoom].TryGetValue(p, out ret))
121
            {
122
               ret = FromPixelToLatLng(p.X, p.Y, zoom);
123
               FromPixelToLatLngCache[zoom].Add(p, ret);
124
 
125
               // for reverse cache
126
               if(!FromLatLngToPixelCache[zoom].ContainsKey(ret))
127
               {
128
                  FromLatLngToPixelCache[zoom].Add(ret, p);
129
               }
130
 
131
               Debug.WriteLine("FromPixelToLatLngCache[" + zoom + "] added " + p + " with " + ret);
132
            }
133
            return ret;
134
         }
135
         else
136
         {
137
            return FromPixelToLatLng(p.X, p.Y, zoom);
138
         }
139
      }
140
 
141
      /// <summary>
142
      /// gets tile coorddinate from pixel coordinates
143
      /// </summary>
144
      /// <param name="p"></param>
145
      /// <returns></returns>
146
      public virtual GPoint FromPixelToTileXY(GPoint p)
147
      {
148
         return new GPoint((long)(p.X / TileSize.Width), (long)(p.Y / TileSize.Height));
149
      }
150
 
151
      /// <summary>
152
      /// gets pixel coordinate from tile coordinate
153
      /// </summary>
154
      /// <param name="p"></param>
155
      /// <returns></returns>
156
      public virtual GPoint FromTileXYToPixel(GPoint p)
157
      {
158
         return new GPoint((p.X * TileSize.Width), (p.Y * TileSize.Height));
159
      }
160
 
161
      /// <summary>
162
      /// min. tile in tiles at custom zoom level
163
      /// </summary>
164
      /// <param name="zoom"></param>
165
      /// <returns></returns>
166
      public abstract GSize GetTileMatrixMinXY(int zoom);
167
 
168
      /// <summary>
169
      /// max. tile in tiles at custom zoom level
170
      /// </summary>
171
      /// <param name="zoom"></param>
172
      /// <returns></returns>
173
      public abstract GSize GetTileMatrixMaxXY(int zoom);
174
 
175
      /// <summary>
176
      /// gets matrix size in tiles
177
      /// </summary>
178
      /// <param name="zoom"></param>
179
      /// <returns></returns>
180
      public virtual GSize GetTileMatrixSizeXY(int zoom)
181
      {
182
         GSize sMin = GetTileMatrixMinXY(zoom);
183
         GSize sMax = GetTileMatrixMaxXY(zoom);
184
 
185
         return new GSize(sMax.Width - sMin.Width + 1, sMax.Height - sMin.Height + 1);
186
      }
187
 
188
      /// <summary>
189
      /// tile matrix size in pixels at custom zoom level
190
      /// </summary>
191
      /// <param name="zoom"></param>
192
      /// <returns></returns>
193
      public long GetTileMatrixItemCount(int zoom)
194
      {
195
         GSize s = GetTileMatrixSizeXY(zoom);
196
         return (s.Width * s.Height);
197
      }
198
 
199
      /// <summary>
200
      /// gets matrix size in pixels
201
      /// </summary>
202
      /// <param name="zoom"></param>
203
      /// <returns></returns>
204
      public virtual GSize GetTileMatrixSizePixel(int zoom)
205
      {
206
         GSize s = GetTileMatrixSizeXY(zoom);
207
         return new GSize(s.Width * TileSize.Width, s.Height * TileSize.Height);
208
      }
209
 
210
      /// <summary>
211
      /// gets all tiles in rect at specific zoom
212
      /// </summary>
213
      public List<GPoint> GetAreaTileList(RectLatLng rect, int zoom, int padding)
214
      {
215
         List<GPoint> ret = new List<GPoint>();
216
 
217
         GPoint topLeft = FromPixelToTileXY(FromLatLngToPixel(rect.LocationTopLeft, zoom));
218
         GPoint rightBottom = FromPixelToTileXY(FromLatLngToPixel(rect.LocationRightBottom, zoom));
219
 
220
         for(long x = (topLeft.X - padding); x <= (rightBottom.X + padding); x++)
221
         {
222
            for(long y = (topLeft.Y - padding); y <= (rightBottom.Y + padding); y++)
223
            {
224
               GPoint p = new GPoint(x, y);
225
               if(!ret.Contains(p) && p.X >= 0 && p.Y >= 0)
226
               {
227
                  ret.Add(p);
228
               }
229
            }
230
         }
231
 
232
         return ret;
233
      }
234
 
235
      /// <summary>
236
      /// The ground resolution indicates the distance (in meters) on the ground that’s represented by a single pixel in the map.
237
      /// For example, at a ground resolution of 10 meters/pixel, each pixel represents a ground distance of 10 meters.
238
      /// </summary>
239
      /// <param name="zoom"></param>
240
      /// <param name="latitude"></param>
241
      /// <returns></returns>
242
      public virtual double GetGroundResolution(int zoom, double latitude)
243
      {
244
         return (Math.Cos(latitude * (Math.PI / 180)) * 2 * Math.PI * Axis) / GetTileMatrixSizePixel(zoom).Width;
245
      }
246
 
247
      /// <summary>
248
      /// gets boundaries
249
      /// </summary>
250
      public virtual RectLatLng Bounds
251
      {
252
         get
253
         {
254
            return RectLatLng.FromLTRB(-180, 90, 180, -90);
255
         }
256
      }
257
 
258
      #region -- math functions --
259
 
260
      /// <summary>
261
      /// PI
262
      /// </summary>
263
      protected static readonly double PI = Math.PI;
264
 
265
      /// <summary>
266
      /// Half of PI
267
      /// </summary>
268
      protected static readonly double HALF_PI = (PI * 0.5);
269
 
270
      /// <summary>
271
      /// PI * 2
272
      /// </summary>
273
      protected static readonly double TWO_PI = (PI * 2.0);
274
 
275
      /// <summary>
276
      /// EPSLoN
277
      /// </summary>
278
      protected static readonly double EPSLoN = 1.0e-10;
279
 
280
      /// <summary>
281
      /// MAX_VAL
282
      /// </summary>
283
      protected const double MAX_VAL = 4;
284
 
285
      /// <summary>
286
      /// MAXLONG
287
      /// </summary>
288
      protected static readonly double MAXLONG = 2147483647;
289
 
290
      /// <summary>
291
      /// DBLLONG
292
      /// </summary>
293
      protected static readonly double DBLLONG = 4.61168601e18;
294
 
295
      static readonly double R2D = 180 / Math.PI;
296
      static readonly double D2R = Math.PI / 180;
297
 
298
      public static double DegreesToRadians(double deg)
299
      {
300
         return (D2R * deg);
301
      }
302
 
303
      public static double RadiansToDegrees(double rad)
304
      {
305
         return (R2D * rad);
306
      }
307
 
308
      ///<summary>
309
      /// return the sign of an argument 
310
      ///</summary>
311
      protected static double Sign(double x)
312
      {
313
         if(x < 0.0)
314
            return (-1);
315
         else
316
            return (1);
317
      }
318
 
319
      protected static double AdjustLongitude(double x)
320
      {
321
         long count = 0;
322
         while(true)
323
         {
324
            if(Math.Abs(x) <= PI)
325
               break;
326
            else
327
               if(((long)Math.Abs(x / Math.PI)) < 2)
328
                  x = x - (Sign(x) * TWO_PI);
329
               else
330
                  if(((long)Math.Abs(x / TWO_PI)) < MAXLONG)
331
                  {
332
                     x = x - (((long)(x / TWO_PI)) * TWO_PI);
333
                  }
334
                  else
335
                     if(((long)Math.Abs(x / (MAXLONG * TWO_PI))) < MAXLONG)
336
                     {
337
                        x = x - (((long)(x / (MAXLONG * TWO_PI))) * (TWO_PI * MAXLONG));
338
                     }
339
                     else
340
                        if(((long)Math.Abs(x / (DBLLONG * TWO_PI))) < MAXLONG)
341
                        {
342
                           x = x - (((long)(x / (DBLLONG * TWO_PI))) * (TWO_PI * DBLLONG));
343
                        }
344
                        else
345
                           x = x - (Sign(x) * TWO_PI);
346
            count++;
347
            if(count > MAX_VAL)
348
               break;
349
         }
350
         return (x);
351
      }
352
 
353
      /// <summary>
354
      /// calculates the sine and cosine
355
      /// </summary>
356
      protected static void SinCos(double val, out double sin, out double cos)
357
      {
358
         sin = Math.Sin(val);
359
         cos = Math.Cos(val);
360
      }
361
 
362
      /// <summary>
363
      /// computes the constants e0, e1, e2, and e3 which are used
364
      /// in a series for calculating the distance along a meridian.
365
      /// </summary>
366
      /// <param name="x">represents the eccentricity squared</param>
367
      /// <returns></returns>
368
      protected static double e0fn(double x)
369
      {
370
         return (1.0 - 0.25 * x * (1.0 + x / 16.0 * (3.0 + 1.25 * x)));
371
      }
372
 
373
      protected static double e1fn(double x)
374
      {
375
         return (0.375 * x * (1.0 + 0.25 * x * (1.0 + 0.46875 * x)));
376
      }
377
 
378
      protected static double e2fn(double x)
379
      {
380
         return (0.05859375 * x * x * (1.0 + 0.75 * x));
381
      }
382
 
383
      protected static double e3fn(double x)
384
      {
385
         return (x * x * x * (35.0 / 3072.0));
386
      }
387
 
388
      /// <summary>
389
      /// computes the value of M which is the distance along a meridian
390
      /// from the Equator to latitude phi.
391
      /// </summary>
392
      protected static double mlfn(double e0, double e1, double e2, double e3, double phi)
393
      {
394
         return (e0 * phi - e1 * Math.Sin(2.0 * phi) + e2 * Math.Sin(4.0 * phi) - e3 * Math.Sin(6.0 * phi));
395
      }
396
 
397
      /// <summary>
398
      /// calculates UTM zone number
399
      /// </summary>
400
      /// <param name="lon">Longitude in degrees</param>
401
      /// <returns></returns>
402
      protected static long GetUTMzone(double lon)
403
      {
404
         return ((long)(((lon + 180.0) / 6.0) + 1.0));
405
      }
406
 
407
      /// <summary>
408
      /// Clips a number to the specified minimum and maximum values.
409
      /// </summary>
410
      /// <param name="n">The number to clip.</param>
411
      /// <param name="minValue">Minimum allowable value.</param>
412
      /// <param name="maxValue">Maximum allowable value.</param>
413
      /// <returns>The clipped value.</returns>
414
      protected static double Clip(double n, double minValue, double maxValue)
415
      {
416
         return Math.Min(Math.Max(n, minValue), maxValue);
417
      }
418
 
419
      /// <summary>
420
      /// distance (in km) between two points specified by latitude/longitude
421
      /// The Haversine formula, http://www.movable-type.co.uk/scripts/latlong.html
422
      /// </summary>
423
      /// <param name="p1"></param>
424
      /// <param name="p2"></param>
425
      /// <returns></returns>
426
      public double GetDistance(PointLatLng p1, PointLatLng p2)
427
      {
428
         double dLat1InRad = p1.Lat * (Math.PI / 180);
429
         double dLong1InRad = p1.Lng * (Math.PI / 180);
430
         double dLat2InRad = p2.Lat * (Math.PI / 180);
431
         double dLong2InRad = p2.Lng * (Math.PI / 180);
432
         double dLongitude = dLong2InRad - dLong1InRad;
433
         double dLatitude = dLat2InRad - dLat1InRad;
434
         double a = Math.Pow(Math.Sin(dLatitude / 2), 2) + Math.Cos(dLat1InRad) * Math.Cos(dLat2InRad) * Math.Pow(Math.Sin(dLongitude / 2), 2);
435
         double c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
436
         double dDistance = (Axis / 1000.0) * c;
437
         return dDistance;
438
      }
439
 
2407 - 440
 
2287 - 441
      public double GetDistanceInPixels(GPoint point1, GPoint point2)
442
      {
443
         double a = (double)(point2.X - point1.X);
444
         double b = (double)(point2.Y - point1.Y);
445
 
446
         return Math.Sqrt(a * a + b * b);
447
      }
448
 
2407 - 449
        /// <summary>
450
        /// point from distance (in m) to point specified by latitude/longitude
451
        /// without bearing (bearing = 0°)
452
        /// for calculating radius of circle around a point
453
        /// The Haversine formula, http://www.movable-type.co.uk/scripts/latlong.html
454
        /// </summary>
455
        /// <param name="p1"></param>
456
        /// <param name="dDist"></param>
457
        /// <returns></returns>
458
        public static PointLatLng GetPointFromDistance(PointLatLng p1, double dDist)
459
      {
460
            double R = 6378137;
461
            double dBearing = 0;
462
            double dLat1InRad = p1.Lat * (Math.PI / 180);
463
            double dLong1InRad = p1.Lng * (Math.PI / 180);
464
 
465
            double dLat2InRad = Math.Asin(Math.Sin(dLat1InRad) * Math.Cos(dDist / R) + Math.Cos(dLat1InRad) * Math.Sin(dDist / R));
466
            double dLong2InRad = dLong1InRad + Math.Atan2(Math.Sin(dBearing) * Math.Sin(dDist / R) * Math.Cos(dLat1InRad), Math.Cos(dDist / R) - Math.Sin(dLat1InRad) * Math.Sin(dLat2InRad));
467
 
468
            double dLatitude = dLat2InRad / (Math.PI / 180);
469
            double dLongitude = dLong2InRad / (Math.PI / 180);
470
 
471
            return new PointLatLng(dLatitude, dLongitude);
472
      }
2287 - 473
      /// <summary>
474
      /// Accepts two coordinates in degrees.
475
      /// </summary>
476
      /// <returns>A double value in degrees. From 0 to 360.</returns>
477
      public double GetBearing(PointLatLng p1, PointLatLng p2)
478
      {
479
         var latitude1 = DegreesToRadians(p1.Lat);
480
         var latitude2 = DegreesToRadians(p2.Lat);
481
         var longitudeDifference = DegreesToRadians(p2.Lng - p1.Lng);
482
 
483
         var y = Math.Sin(longitudeDifference) * Math.Cos(latitude2);
484
         var x = Math.Cos(latitude1) * Math.Sin(latitude2) - Math.Sin(latitude1) * Math.Cos(latitude2) * Math.Cos(longitudeDifference);
485
 
486
         return (RadiansToDegrees(Math.Atan2(y, x)) + 360) % 360;
487
      }
488
 
489
      /// <summary>
490
      /// Conversion from cartesian earth-sentered coordinates to geodetic coordinates in the given datum
491
      /// </summary>
492
      /// <param name="Lat"></param>
493
      /// <param name="Lon"></param>
494
      /// <param name="Height">Height above ellipsoid [m]</param>
495
      /// <param name="X"></param>
496
      /// <param name="Y"></param>
497
      /// <param name="Z"></param>
498
      public void FromGeodeticToCartesian(double Lat, double Lng, double Height, out double X, out double Y, out double Z)
499
      {
500
         Lat = (Math.PI / 180) * Lat;
501
         Lng = (Math.PI / 180) * Lng;
502
 
503
         double B = Axis * (1.0 - Flattening);
504
         double ee = 1.0 - (B / Axis) * (B / Axis);
505
         double N = (Axis / Math.Sqrt(1.0 - ee * Math.Sin(Lat) * Math.Sin(Lat)));
506
 
507
         X = (N + Height) * Math.Cos(Lat) * Math.Cos(Lng);
508
         Y = (N + Height) * Math.Cos(Lat) * Math.Sin(Lng);
509
         Z = (N * (B / Axis) * (B / Axis) + Height) * Math.Sin(Lat);
510
      }
511
 
512
      /// <summary>
513
      /// Conversion from cartesian earth-sentered coordinates to geodetic coordinates in the given datum
514
      /// </summary>
515
      /// <param name="X"></param>
516
      /// <param name="Y"></param>
517
      /// <param name="Z"></param>
518
      /// <param name="Lat"></param>
519
      /// <param name="Lon"></param>
520
      public void FromCartesianTGeodetic(double X, double Y, double Z, out double Lat, out double Lng)
521
      {
522
         double E = Flattening * (2.0 - Flattening);
523
         Lng = Math.Atan2(Y, X);
524
 
525
         double P = Math.Sqrt(X * X + Y * Y);
526
         double Theta = Math.Atan2(Z, (P * (1.0 - Flattening)));
527
         double st = Math.Sin(Theta);
528
         double ct = Math.Cos(Theta);
529
         Lat = Math.Atan2(Z + E / (1.0 - Flattening) * Axis * st * st * st, P - E * Axis * ct * ct * ct);
530
 
531
         Lat /= (Math.PI / 180);
532
         Lng /= (Math.PI / 180);
533
      }
534
 
535
      #endregion
536
   }
537
}