Subversion Repositories Projects

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2498 - 1

2
namespace GMap.NET.Projections
3
{
4
    using System;
5
    using System.Collections.Generic;
6
 
7
    /// <summary>
8
    /// PROJCS["SWEREF99 TM",GEOGCS["SWEREF99",DATUM["SWEREF99",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6619"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4619"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],
9
    /// PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],AUTHORITY["EPSG","3006"],AXIS["y",EAST],AXIS["x",NORTH]]
10
    /// </summary>
11
    public class SWEREF99_TMProjection : PureProjection
12
    {
13
        public static readonly SWEREF99_TMProjection Instance = new SWEREF99_TMProjection();
14
 
15
        static readonly double MinLatitude = 54.96;
16
        static readonly double MaxLatitude = 69.07;
17
        static readonly double MinLongitude = 10.0;
18
        static readonly double MaxLongitude = 24.5;
19
 
20
        static readonly double orignX = -1200000;
21
        static readonly double orignY = 8500000;
22
 
23
        static readonly double scaleFactor = 0.9996;                    // scale factor                         
24
        static readonly double centralMeridian = DegreesToRadians(15);// Center longitude (projection center) 
25
        static readonly double latOrigin = 0.0;                    // center latitude                   
26
        static readonly double falseNorthing = 0.0;           // y offset in meters                     
27
        static readonly double falseEasting = 500000.0;        // x offset in meters                    
28
        static readonly double semiMajor = 6378137.0;           // major axis
29
        static readonly double semiMinor = 6356752.3141403561; // minor axis
30
        static readonly double semiMinor2 = 6356752.3142451793;     // minor axis
31
        static readonly double metersPerUnit = 1.0;
32
        static readonly double COS_67P5 = 0.38268343236508977; // cosine of 67.5 degrees
33
        static readonly double AD_C = 1.0026000;               // Toms region 1 constant
34
 
35
        public override RectLatLng Bounds
36
        {
37
            get
38
            {
39
                return RectLatLng.FromLTRB(MinLongitude, MaxLatitude, MaxLongitude, MinLatitude);
40
            }
41
        }
42
 
43
        GSize tileSize = new GSize(256, 256);
44
        public override GSize TileSize
45
        {
46
            get
47
            {
48
                return tileSize;
49
            }
50
        }
51
 
52
        public override double Axis
53
        {
54
            get
55
            {
56
                return 6378137;
57
            }
58
        }
59
 
60
        public override double Flattening
61
        {
62
            get
63
            {
64
                return (1.0 / 298.257222101);
65
            }
66
        }
67
 
68
        public override GPoint FromLatLngToPixel(double lat, double lng, int zoom)
69
        {
70
            GPoint ret = GPoint.Empty;
71
 
72
            lat = Clip(lat, MinLatitude, MaxLatitude);
73
            lng = Clip(lng, MinLongitude, MaxLongitude);
74
 
75
            double[] lks = new double[] { lng, lat };
76
            lks = DTM10(lks);
77
            lks = MTD10(lks);
78
            lks = DTM00(lks);
79
 
80
            double res = GetTileMatrixResolution(zoom);
81
            return LksToPixel(lks, res);
82
        }
83
 
84
        static GPoint LksToPixel(double[] lks, double res)
85
        {
86
            return new GPoint((long)Math.Floor((lks[0] - orignX) / res), (long)Math.Floor((orignY - lks[1]) / res));
87
        }
88
 
89
        public override PointLatLng FromPixelToLatLng(long x, long y, int zoom)
90
        {
91
            PointLatLng ret = PointLatLng.Empty;
92
 
93
            double res = GetTileMatrixResolution(zoom);
94
 
95
            double[] lks = new double[] { (x * res) + orignX, orignY - (y * res) };
96
            lks = MTD11(lks);
97
            lks = DTM10(lks);
98
            lks = MTD10(lks);
99
 
100
            ret.Lat = Clip(lks[1], MinLatitude, MaxLatitude);
101
            ret.Lng = Clip(lks[0], MinLongitude, MaxLongitude);
102
 
103
            return ret;
104
        }
105
 
106
        double[] DTM10(double[] lonlat)
107
        {
108
            // Eccentricity squared : (a^2 - b^2)/a^2
109
            double es = 1.0 - (semiMinor2 * semiMinor2) / (semiMajor * semiMajor); // e^2
110
 
111
            // Second eccentricity squared : (a^2 - b^2)/b^2
112
            double ses = (Math.Pow(semiMajor, 2) - Math.Pow(semiMinor2, 2)) / Math.Pow(semiMinor2, 2);
113
 
114
            double ba = semiMinor2 / semiMajor;
115
            double ab = semiMajor / semiMinor2;
116
 
117
            double lon = DegreesToRadians(lonlat[0]);
118
            double lat = DegreesToRadians(lonlat[1]);
119
            double h = lonlat.Length < 3 ? 0 : lonlat[2].Equals(Double.NaN) ? 0 : lonlat[2];
120
            double v = semiMajor / Math.Sqrt(1 - es * Math.Pow(Math.Sin(lat), 2));
121
            double x = (v + h) * Math.Cos(lat) * Math.Cos(lon);
122
            double y = (v + h) * Math.Cos(lat) * Math.Sin(lon);
123
            double z = ((1 - es) * v + h) * Math.Sin(lat);
124
            return new double[] { x, y, z, };
125
        }
126
 
127
        double[] MTD10(double[] pnt)
128
        {
129
            // Eccentricity squared : (a^2 - b^2)/a^2
130
            double es = 1.0 - (semiMinor * semiMinor) / (semiMajor * semiMajor); // e^2
131
 
132
            // Second eccentricity squared : (a^2 - b^2)/b^2
133
            double ses = (Math.Pow(semiMajor, 2) - Math.Pow(semiMinor, 2)) / Math.Pow(semiMinor, 2);
134
 
135
            double ba = semiMinor / semiMajor;
136
            double ab = semiMajor / semiMinor;
137
 
138
            bool AtPole = false; // is location in polar region
139
            double Z = pnt.Length < 3 ? 0 : pnt[2].Equals(Double.NaN) ? 0 : pnt[2];
140
 
141
            double lon = 0;
142
            double lat = 0;
143
            double Height = 0;
144
            if (pnt[0] != 0.0)
145
            {
146
                lon = Math.Atan2(pnt[1], pnt[0]);
147
            }
148
            else
149
            {
150
                if (pnt[1] > 0)
151
                {
152
                    lon = Math.PI / 2;
153
                }
154
                else
155
                   if (pnt[1] < 0)
156
                {
157
                    lon = -Math.PI * 0.5;
158
                }
159
                else
160
                {
161
                    AtPole = true;
162
                    lon = 0.0;
163
                    if (Z > 0.0) // north pole
164
                    {
165
                        lat = Math.PI * 0.5;
166
                    }
167
                    else
168
                       if (Z < 0.0) // south pole
169
                    {
170
                        lat = -Math.PI * 0.5;
171
                    }
172
                    else // center of earth
173
                    {
174
                        return new double[] { RadiansToDegrees(lon), RadiansToDegrees(Math.PI * 0.5), -semiMinor, };
175
                    }
176
                }
177
            }
178
            double W2 = pnt[0] * pnt[0] + pnt[1] * pnt[1]; // Square of distance from Z axis
179
            double W = Math.Sqrt(W2); // distance from Z axis
180
            double T0 = Z * AD_C; // initial estimate of vertical component
181
            double S0 = Math.Sqrt(T0 * T0 + W2); // initial estimate of horizontal component
182
            double Sin_B0 = T0 / S0; // sin(B0), B0 is estimate of Bowring aux variable
183
            double Cos_B0 = W / S0; // cos(B0)
184
            double Sin3_B0 = Math.Pow(Sin_B0, 3);
185
            double T1 = Z + semiMinor * ses * Sin3_B0; // corrected estimate of vertical component
186
            double Sum = W - semiMajor * es * Cos_B0 * Cos_B0 * Cos_B0; // numerator of cos(phi1)
187
            double S1 = Math.Sqrt(T1 * T1 + Sum * Sum); // corrected estimate of horizontal component
188
            double Sin_p1 = T1 / S1; // sin(phi1), phi1 is estimated latitude
189
            double Cos_p1 = Sum / S1; // cos(phi1)
190
            double Rn = semiMajor / Math.Sqrt(1.0 - es * Sin_p1 * Sin_p1); // Earth radius at location
191
            if (Cos_p1 >= COS_67P5)
192
            {
193
                Height = W / Cos_p1 - Rn;
194
            }
195
            else
196
               if (Cos_p1 <= -COS_67P5)
197
            {
198
                Height = W / -Cos_p1 - Rn;
199
            }
200
            else
201
            {
202
                Height = Z / Sin_p1 + Rn * (es - 1.0);
203
            }
204
 
205
            if (!AtPole)
206
            {
207
                lat = Math.Atan(Sin_p1 / Cos_p1);
208
            }
209
            return new double[] { RadiansToDegrees(lon), RadiansToDegrees(lat), Height, };
210
        }
211
 
212
        double[] DTM00(double[] lonlat)
213
        {
214
            double e0, e1, e2, e3;  // eccentricity constants           
215
            double e, es, esp;      // eccentricity constants           
216
            double ml0;              // small value m                   
217
 
218
            es = 1.0 - Math.Pow(semiMinor / semiMajor, 2);
219
            e = Math.Sqrt(es);
220
            e0 = e0fn(es);
221
            e1 = e1fn(es);
222
            e2 = e2fn(es);
223
            e3 = e3fn(es);
224
            ml0 = semiMajor * mlfn(e0, e1, e2, e3, latOrigin);
225
            esp = es / (1.0 - es);
226
 
227
            // ...              
228
 
229
            double lon = DegreesToRadians(lonlat[0]);
230
            double lat = DegreesToRadians(lonlat[1]);
231
 
232
            double delta_lon = 0.0;  // Delta longitude (Given longitude - center)
233
            double sin_phi, cos_phi; // sin and cos value                               
234
            double al, als;         // temporary values                         
235
            double c, t, tq;           // temporary values                              
236
            double con, n, ml;      // cone constant, small m                   
237
 
238
            delta_lon = AdjustLongitude(lon - centralMeridian);
239
            SinCos(lat, out sin_phi, out cos_phi);
240
 
241
            al = cos_phi * delta_lon;
242
            als = Math.Pow(al, 2);
243
            c = esp * Math.Pow(cos_phi, 2);
244
            tq = Math.Tan(lat);
245
            t = Math.Pow(tq, 2);
246
            con = 1.0 - es * Math.Pow(sin_phi, 2);
247
            n = semiMajor / Math.Sqrt(con);
248
            ml = semiMajor * mlfn(e0, e1, e2, e3, lat);
249
 
250
            double x = scaleFactor * n * al * (1.0 + als / 6.0 * (1.0 - t + c + als / 20.0 *
251
                (5.0 - 18.0 * t + Math.Pow(t, 2) + 72.0 * c - 58.0 * esp))) + falseEasting;
252
 
253
            double y = scaleFactor * (ml - ml0 + n * tq * (als * (0.5 + als / 24.0 *
254
                (5.0 - t + 9.0 * c + 4.0 * Math.Pow(c, 2) + als / 30.0 * (61.0 - 58.0 * t
255
                + Math.Pow(t, 2) + 600.0 * c - 330.0 * esp))))) + falseNorthing;
256
 
257
            if (lonlat.Length < 3)
258
                return new double[] { x / metersPerUnit, y / metersPerUnit };
259
            else
260
                return new double[] { x / metersPerUnit, y / metersPerUnit, lonlat[2] };
261
        }
262
 
263
        double[] DTM01(double[] lonlat)
264
        {
265
            // Eccentricity squared : (a^2 - b^2)/a^2
266
            double es = 1.0 - (semiMinor * semiMinor) / (semiMajor * semiMajor);
267
 
268
            // Second eccentricity squared : (a^2 - b^2)/b^2
269
            double ses = (Math.Pow(semiMajor, 2) - Math.Pow(semiMinor, 2)) / Math.Pow(semiMinor, 2);
270
 
271
            double ba = semiMinor / semiMajor;
272
            double ab = semiMajor / semiMinor;
273
 
274
            double lon = DegreesToRadians(lonlat[0]);
275
            double lat = DegreesToRadians(lonlat[1]);
276
            double h = lonlat.Length < 3 ? 0 : lonlat[2].Equals(Double.NaN) ? 0 : lonlat[2];
277
            double v = semiMajor / Math.Sqrt(1 - es * Math.Pow(Math.Sin(lat), 2));
278
            double x = (v + h) * Math.Cos(lat) * Math.Cos(lon);
279
            double y = (v + h) * Math.Cos(lat) * Math.Sin(lon);
280
            double z = ((1 - es) * v + h) * Math.Sin(lat);
281
            return new double[] { x, y, z, };
282
        }
283
 
284
        double[] MTD01(double[] pnt)
285
        {
286
            // Eccentricity squared : (a^2 - b^2)/a^2
287
            double es = 1.0 - (semiMinor2 * semiMinor2) / (semiMajor * semiMajor);
288
 
289
            // Second eccentricity squared : (a^2 - b^2)/b^2
290
            double ses = (Math.Pow(semiMajor, 2) - Math.Pow(semiMinor2, 2)) / Math.Pow(semiMinor2, 2);
291
 
292
            double ba = semiMinor2 / semiMajor;
293
            double ab = semiMajor / semiMinor2;
294
 
295
            bool At_Pole = false; // is location in polar region
296
            double Z = pnt.Length < 3 ? 0 : pnt[2].Equals(Double.NaN) ? 0 : pnt[2];
297
 
298
            double lon = 0;
299
            double lat = 0;
300
            double Height = 0;
301
            if (pnt[0] != 0.0)
302
            {
303
                lon = Math.Atan2(pnt[1], pnt[0]);
304
            }
305
            else
306
            {
307
                if (pnt[1] > 0)
308
                {
309
                    lon = Math.PI / 2;
310
                }
311
                else
312
                   if (pnt[1] < 0)
313
                {
314
                    lon = -Math.PI * 0.5;
315
                }
316
                else
317
                {
318
                    At_Pole = true;
319
                    lon = 0.0;
320
                    if (Z > 0.0) // north pole
321
                    {
322
                        lat = Math.PI * 0.5;
323
                    }
324
                    else
325
                       if (Z < 0.0) // south pole
326
                    {
327
                        lat = -Math.PI * 0.5;
328
                    }
329
                    else // center of earth
330
                    {
331
                        return new double[] { RadiansToDegrees(lon), RadiansToDegrees(Math.PI * 0.5), -semiMinor2, };
332
                    }
333
                }
334
            }
335
 
336
            double W2 = pnt[0] * pnt[0] + pnt[1] * pnt[1]; // Square of distance from Z axis
337
            double W = Math.Sqrt(W2);                      // distance from Z axis
338
            double T0 = Z * AD_C;                // initial estimate of vertical component
339
            double S0 = Math.Sqrt(T0 * T0 + W2); //initial estimate of horizontal component
340
            double Sin_B0 = T0 / S0;             // sin(B0), B0 is estimate of Bowring aux variable
341
            double Cos_B0 = W / S0;              // cos(B0)
342
            double Sin3_B0 = Math.Pow(Sin_B0, 3);
343
            double T1 = Z + semiMinor2 * ses * Sin3_B0; //corrected estimate of vertical component
344
            double Sum = W - semiMajor * es * Cos_B0 * Cos_B0 * Cos_B0; // numerator of cos(phi1)
345
            double S1 = Math.Sqrt(T1 * T1 + Sum * Sum); // corrected estimate of horizontal component
346
            double Sin_p1 = T1 / S1;  // sin(phi1), phi1 is estimated latitude
347
            double Cos_p1 = Sum / S1; // cos(phi1)
348
            double Rn = semiMajor / Math.Sqrt(1.0 - es * Sin_p1 * Sin_p1); // Earth radius at location
349
 
350
            if (Cos_p1 >= COS_67P5)
351
            {
352
                Height = W / Cos_p1 - Rn;
353
            }
354
            else
355
               if (Cos_p1 <= -COS_67P5)
356
            {
357
                Height = W / -Cos_p1 - Rn;
358
            }
359
            else
360
            {
361
                Height = Z / Sin_p1 + Rn * (es - 1.0);
362
            }
363
 
364
            if (!At_Pole)
365
            {
366
                lat = Math.Atan(Sin_p1 / Cos_p1);
367
            }
368
            return new double[] { RadiansToDegrees(lon), RadiansToDegrees(lat), Height, };
369
        }
370
 
371
        double[] MTD11(double[] p)
372
        {
373
            double e0, e1, e2, e3;  // eccentricity constants           
374
            double e, es, esp;      // eccentricity constants           
375
            double ml0;         // small value m
376
 
377
            es = 1.0 - Math.Pow(semiMinor / semiMajor, 2);
378
            e = Math.Sqrt(es);
379
            e0 = e0fn(es);
380
            e1 = e1fn(es);
381
            e2 = e2fn(es);
382
            e3 = e3fn(es);
383
            ml0 = semiMajor * mlfn(e0, e1, e2, e3, latOrigin);
384
            esp = es / (1.0 - es);
385
 
386
            // ...
387
 
388
            double con, phi;
389
            double delta_phi;
390
            long i;
391
            double sin_phi, cos_phi, tan_phi;
392
            double c, cs, t, ts, n, r, d, ds;
393
            long max_iter = 6;
394
 
395
            double x = p[0] * metersPerUnit - falseEasting;
396
            double y = p[1] * metersPerUnit - falseNorthing;
397
 
398
            con = (ml0 + y / scaleFactor) / semiMajor;
399
            phi = con;
400
            for (i = 0; ; i++)
401
            {
402
                delta_phi = ((con + e1 * Math.Sin(2.0 * phi) - e2 * Math.Sin(4.0 * phi) + e3 * Math.Sin(6.0 * phi)) / e0) - phi;
403
                phi += delta_phi;
404
 
405
                if (Math.Abs(delta_phi) <= EPSLoN)
406
                    break;
407
 
408
                if (i >= max_iter)
409
                    throw new ArgumentException("Latitude failed to converge");
410
            }
411
 
412
            if (Math.Abs(phi) < HALF_PI)
413
            {
414
                SinCos(phi, out sin_phi, out cos_phi);
415
                tan_phi = Math.Tan(phi);
416
                c = esp * Math.Pow(cos_phi, 2);
417
                cs = Math.Pow(c, 2);
418
                t = Math.Pow(tan_phi, 2);
419
                ts = Math.Pow(t, 2);
420
                con = 1.0 - es * Math.Pow(sin_phi, 2);
421
                n = semiMajor / Math.Sqrt(con);
422
                r = n * (1.0 - es) / con;
423
                d = x / (n * scaleFactor);
424
                ds = Math.Pow(d, 2);
425
 
426
                double lat = phi - (n * tan_phi * ds / r) * (0.5 - ds / 24.0 * (5.0 + 3.0 * t +
427
                    10.0 * c - 4.0 * cs - 9.0 * esp - ds / 30.0 * (61.0 + 90.0 * t +
428
                    298.0 * c + 45.0 * ts - 252.0 * esp - 3.0 * cs)));
429
 
430
                double lon = AdjustLongitude(centralMeridian + (d * (1.0 - ds / 6.0 * (1.0 + 2.0 * t +
431
                    c - ds / 20.0 * (5.0 - 2.0 * c + 28.0 * t - 3.0 * cs + 8.0 * esp +
432
                    24.0 * ts))) / cos_phi));
433
 
434
                if (p.Length < 3)
435
                    return new double[] { RadiansToDegrees(lon), RadiansToDegrees(lat) };
436
                else
437
                    return new double[] { RadiansToDegrees(lon), RadiansToDegrees(lat), p[2] };
438
            }
439
            else
440
            {
441
                if (p.Length < 3)
442
                    return new double[] { RadiansToDegrees(HALF_PI * Sign(y)), RadiansToDegrees(centralMeridian) };
443
                else
444
                    return new double[] { RadiansToDegrees(HALF_PI * Sign(y)), RadiansToDegrees(centralMeridian), p[2] };
445
            }
446
        }
447
 
448
        #region -- levels info --
449
        //      "defaultLayer" : "topowebbwmts",
450
        //      "extent" : {
451
        //              "left" : -1200000,
452
        //              "bottom" : 4700000,
453
        //              "right" : 2600000,
454
        //              "top" : 8500000
455
        //      },
456
        //      "projection" : "EPSG:3006",
457
        //      "units" : "m",
458
        //      "allOverlays" : true,
459
        //      "resolutions" : [4096.0, 2048.0, 1024.0, 512.0, 256.0, 128.0, 64.0, 32.0, 16.0, 8.0, 4.0, 2.0, 1.0, 0.5, 0.25, 0.15, 0.1, 0.05, 0.01],
460
        //      "initPosition" : {
461
        //              "n" : 6607899,
462
        //              "e" : 564931,
463
        //              "zoom" : 2
464
        //      },
465
        #endregion
466
 
467
        static double[] resolutions = new double[] { 4096.0, 2048.0, 1024.0, 512.0, 256.0, 128.0, 64.0, 32.0, 16.0, 8.0, 4.0, 2.0, 1.0, 0.5, 0.25, 0.15, 0.1, 0.05, 0.01 };
468
        public static double GetTileMatrixResolution(int zoom)
469
        {
470
            double ret = 0;
471
 
472
            if (zoom < resolutions.Length)
473
            {
474
                ret = resolutions[zoom];
475
            }
476
 
477
            return ret;
478
        }
479
 
480
        public override double GetGroundResolution(int zoom, double latitude)
481
        {
482
            return GetTileMatrixResolution(zoom);
483
        }
484
 
485
        Dictionary<int, GSize> extentMatrixMin;
486
        Dictionary<int, GSize> extentMatrixMax;
487
 
488
        public override GSize GetTileMatrixMinXY(int zoom)
489
        {
490
            if (extentMatrixMin == null)
491
            {
492
                GenerateExtents();
493
            }
494
            return extentMatrixMin[zoom];
495
        }
496
 
497
        public override GSize GetTileMatrixMaxXY(int zoom)
498
        {
499
            if (extentMatrixMax == null)
500
            {
501
                GenerateExtents();
502
            }
503
            return extentMatrixMax[zoom];
504
        }
505
 
506
        void GenerateExtents()
507
        {
508
            extentMatrixMin = new Dictionary<int, GSize>();
509
            extentMatrixMax = new Dictionary<int, GSize>();
510
 
511
            for (int i = 0; i <= resolutions.Length; i++)
512
            {
513
                double res = GetTileMatrixResolution(i);
514
 
515
                extentMatrixMin.Add(i, new GSize(FromPixelToTileXY(FromLatLngToPixel(Bounds.LocationTopLeft, i))));
516
                extentMatrixMax.Add(i, new GSize(FromPixelToTileXY(FromLatLngToPixel(Bounds.LocationRightBottom, i))));
517
            }
518
        }
519
    }
520
}