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
 
6
   /// <summary>
7
   /// GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]
8
   /// PROJCS["Mapy.cz",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",134400000],PARAMETER["false_northing",-41600000],UNIT["1/32meter",0.03125]]
9
   /// </summary>
10
   public class MapyCZProjection : PureProjection
11
   {
12
      public static readonly MapyCZProjection Instance = new MapyCZProjection();
13
 
14
      static readonly double MinLatitude = 26;
15
      static readonly double MaxLatitude = 76;
16
      static readonly double MinLongitude = -26;
17
      static readonly double MaxLongitude = 38;
18
 
19
      #region -- Common --
20
      static int getLCM(int zone)
21
      {
22
         if((zone < 1) || (zone > 60))
23
         {
24
            throw new Exception("MapyCZProjection: UTM Zone number is not between 1 and 60.");
25
         }
26
         else
27
         {
28
            return ((zone * 6) - 183);
29
         }
30
      }
31
 
32
      static double roundoff(double xx, double yy)
33
      {
34
         var x = xx;
35
         var y = yy;
36
         x = Math.Round(x * Math.Pow(10, y)) / Math.Pow(10, y);
37
         return x;
38
      }
39
 
40
      static readonly double UTMSIZE = 2;
41
      static readonly double UNITS = 1;
42
 
43
      #endregion
44
 
45
      #region -- WGSToMapyCZ --
46
 
47
      public long[] WGSToPP(double la, double lo)
48
      {
49
         var utmEE = wgsToUTM(DegreesToRadians(la), DegreesToRadians(lo), 33);
50
         var pp = utmEEToPP(utmEE[0], utmEE[1]);
51
         return pp;
52
      }
53
 
54
      static long[] utmEEToPP(double east, double north)
55
      {
56
         var x = (Math.Round(east) - (-3700000.0)) * Math.Pow(2, 5);
57
         var y = (Math.Round(north) - (1300000.0)) * Math.Pow(2, 5);
58
 
59
         return new long[] { (long)x, (long)y };
60
      }
61
 
62
      double[] wgsToUTM(double la, double lo, int zone)
63
      {
64
         var latrad = la;
65
         var lonrad = lo;
66
         var latddd = RadiansToDegrees(la);
67
         var londdd = RadiansToDegrees(lo);
68
 
69
         var k = 0.9996f;
70
         var a = Axis;
71
         var f = Flattening;
72
         var b = a * (1.0 - f);
73
         var e2 = (a * a - b * b) / (a * a);
74
         var e = Math.Sqrt(e2);
75
         var ei2 = (a * a - b * b) / (b * b);
76
         var ei = Math.Sqrt(ei2);
77
         var n = (a - b) / (a + b);
78
         var G = a * (1.0 - n) * (1.0 - n * n) * (1.0 + (9 / 4.0) * n * n + (255.0 / 64.0) * Math.Pow(n, 4)) * (PI / 180.0);
79
         var w = londdd - ((double) (zone * 6 - 183));
80
         w = DegreesToRadians(w);
81
         var t = Math.Tan(latrad);
82
         var rho = a * (1.0 - e2) / Math.Pow(1.0 - (e2 * Math.Sin(latrad) * Math.Sin(latrad)), (3 / 2.0));
83
         var nu = a / Math.Sqrt(1.0 - (e2 * Math.Sin(latrad) * Math.Sin(latrad)));
84
         var psi = nu / rho;
85
         var coslat = Math.Cos(latrad);
86
         var sinlat = Math.Sin(latrad);
87
         var A0 = 1 - (e2 / 4.0) - (3 * e2 * e2 / 64.0) - (5 * Math.Pow(e2, 3) / 256.0);
88
         var A2 = (3 / 8.0) * (e2 + (e2 * e2 / 4.0) + (15 * Math.Pow(e2, 3) / 128.0));
89
         var A4 = (15 / 256.0) * (e2 * e2 + (3 * Math.Pow(e2, 3) / 4.0));
90
         var A6 = 35 * Math.Pow(e2, 3) / 3072.0;
91
         var m = a * ((A0 * latrad) - (A2 * Math.Sin(2 * latrad)) + (A4 * Math.Sin(4 * latrad)) - (A6 * Math.Sin(6 * latrad)));
92
         var eterm1 = (w * w / 6.0) * coslat * coslat * (psi - t * t);
93
         var eterm2 = (Math.Pow(w, 4) / 120.0) * Math.Pow(coslat, 4) * (4 * Math.Pow(psi, 3) * (1.0 - 6 * t * t) + psi * psi * (1.0 + 8 * t * t) - psi * 2 * t * t + Math.Pow(t, 4));
94
         var eterm3 = (Math.Pow(w, 6) / 5040.0) * Math.Pow(coslat, 6) * (61.0 - 479 * t * t + 179 * Math.Pow(t, 4) - Math.Pow(t, 6));
95
         var dE = k * nu * w * coslat * (1.0 + eterm1 + eterm2 + eterm3);
96
         var east = 500000.0 + (dE / UNITS);
97
         east = roundoff(east, UTMSIZE);
98
         var nterm1 = (w * w / 2.0) * nu * sinlat * coslat;
99
         var nterm2 = (Math.Pow(w, 4) / 24.0) * nu * sinlat * Math.Pow(coslat, 3) * (4 * psi * psi + psi - t * t);
100
         var nterm3 = (Math.Pow(w, 6) / 720.0) * nu * sinlat * Math.Pow(coslat, 5) * (8 * Math.Pow(psi, 4) * (11.0 - 24 * t * t) - 28 * Math.Pow(psi, 3) * (1.0 - 6 * t * t) + psi * psi * (1.0 - 32 * t * t) - psi * 2 * t * t + Math.Pow(t, 4));
101
         var nterm4 = (Math.Pow(w, 8) / 40320.0) * nu * sinlat * Math.Pow(coslat, 7) * (1385.0 - 3111 * t * t + 543 * Math.Pow(t, 4) - Math.Pow(t, 6));
102
         var dN = k * (m + nterm1 + nterm2 + nterm3 + nterm4);
103
         var north = (0.0 + (dN / UNITS));
104
         north = roundoff(north, UTMSIZE);
105
 
106
         return new double[] { east, north, zone };
107
      }
108
 
109
      #endregion
110
 
111
      #region -- MapyCZToWGS --
112
 
113
      public double[] PPToWGS(double x, double y)
114
      {
115
         var utmEE = ppToUTMEE(x, y);
116
         var ret = utmToWGS(utmEE[0], utmEE[1], 33);
117
         return ret;
118
      }
119
 
120
      double[] ppToUTMEE(double x, double y)
121
      {
122
         var north = y * Math.Pow(2, -5) + 1300000.0;
123
         var east = x * Math.Pow(2, -5) + (-3700000.0);
124
         east = roundoff(east, UTMSIZE);
125
         north = roundoff(north, UTMSIZE);
126
 
127
         return new double[] { east, north };
128
      }
129
 
130
      double[] utmToWGS(double eastIn, double northIn, int zone)
131
      {
132
         var k = 0.9996f;
133
         var a = Axis;
134
         var f = Flattening;
135
         var b = a * (1.0 - f);
136
         var e2 = (a * a - b * b) / (a * a);
137
         var e = Math.Sqrt(e2);
138
         var ei2 = (a * a - b * b) / (b * b);
139
         var ei = Math.Sqrt(ei2);
140
         var n = (a - b) / (a + b);
141
         var G = a * (1.0 - n) * (1.0 - n * n) * (1.0 + (9 / 4.0) * n * n + (255 / 64.0) * Math.Pow(n, 4)) * (PI / 180.0);
142
         var north = (northIn - 0) * UNITS;
143
         var east = (eastIn - 500000.0) * UNITS;
144
         var m = north / k;
145
         var sigma = (m * PI) / (180.0 * G);
146
         var footlat = sigma + ((3 * n / 2.0) - (27 * Math.Pow(n, 3) / 32.0)) * Math.Sin(2 * sigma) + ((21 * n * n / 16.0) - (55 * Math.Pow(n, 4) / 32.0)) * Math.Sin(4 * sigma) + (151 * Math.Pow(n, 3) / 96.0) * Math.Sin(6 * sigma) + (1097 * Math.Pow(n, 4) / 512.0) * Math.Sin(8 * sigma);
147
         var rho = a * (1.0 - e2) / Math.Pow(1.0 - (e2 * Math.Sin(footlat) * Math.Sin(footlat)), (3 / 2.0));
148
         var nu = a / Math.Sqrt(1.0 - (e2 * Math.Sin(footlat) * Math.Sin(footlat)));
149
         var psi = nu / rho;
150
         var t = Math.Tan(footlat);
151
         var x = east / (k * nu);
152
         var laterm1 = (t / (k * rho)) * (east * x / 2.0);
153
         var laterm2 = (t / (k * rho)) * (east * Math.Pow(x, 3) / 24.0) * (-4 * psi * psi + 9 * psi * (1 - t * t) + 12 * t * t);
154
         var laterm3 = (t / (k * rho)) * (east * Math.Pow(x, 5) / 720.0) * (8 * Math.Pow(psi, 4) * (11 - 24 * t * t) - 12 * Math.Pow(psi, 3) * (21.0 - 71 * t * t) + 15 * psi * psi * (15.0 - 98 * t * t + 15 * Math.Pow(t, 4)) + 180 * psi * (5 * t * t - 3 * Math.Pow(t, 4)) + 360 * Math.Pow(t, 4));
155
         var laterm4 = (t / (k * rho)) * (east * Math.Pow(x, 7) / 40320.0) * (1385.0 + 3633 * t * t + 4095 * Math.Pow(t, 4) + 1575 * Math.Pow(t, 6));
156
         var latrad = footlat - laterm1 + laterm2 - laterm3 + laterm4;
157
         var lat = RadiansToDegrees(latrad);
158
         var seclat = 1 / Math.Cos(footlat);
159
         var loterm1 = x * seclat;
160
         var loterm2 = (Math.Pow(x, 3) / 6.0) * seclat * (psi + 2 * t * t);
161
         var loterm3 = (Math.Pow(x, 5) / 120.0) * seclat * (-4 * Math.Pow(psi, 3) * (1 - 6 * t * t) + psi * psi * (9 - 68 * t * t) + 72 * psi * t * t + 24 * Math.Pow(t, 4));
162
         var loterm4 = (Math.Pow(x, 7) / 5040.0) * seclat * (61.0 + 662 * t * t + 1320 * Math.Pow(t, 4) + 720 * Math.Pow(t, 6));
163
         var w = loterm1 - loterm2 + loterm3 - loterm4;
164
         var longrad = DegreesToRadians(getLCM(zone)) + w;
165
         var lon = RadiansToDegrees(longrad);
166
 
167
         return new double[] { lat, lon, latrad, longrad };
168
      }
169
 
170
      #endregion
171
 
172
      public override RectLatLng Bounds
173
      {
174
         get
175
         {
176
            return RectLatLng.FromLTRB(MinLongitude, MaxLatitude, MaxLongitude, MinLatitude);
177
         }
178
      }
179
 
180
      public override GSize TileSize
181
      {
182
         get
183
         {
184
            return new GSize(256, 256);
185
         }
186
      }
187
 
188
      public override double Axis
189
      {
190
         get
191
         {
192
            return 6378137;
193
         }
194
      }
195
 
196
      public override double Flattening
197
      {
198
         get
199
         {
200
            return (1.0 / 298.257223563);
201
         }
202
      }
203
 
204
      public override GPoint FromLatLngToPixel(double lat, double lng, int zoom)
205
      {
206
         GPoint ret = GPoint.Empty;
207
 
208
         lat = Clip(lat, MinLatitude, MaxLatitude);
209
         lng = Clip(lng, MinLongitude, MaxLongitude);
210
 
211
         var size = GetTileMatrixSizePixel(zoom);
212
         {
213
            var l = WGSToPP(lat, lng);
214
            ret.X = (long)l[0] >> (20 - zoom);
215
            ret.Y = size.Height - ((long)l[1] >> (20 - zoom));
216
         }
217
         return ret;
218
      }
219
 
220
      public override PointLatLng FromPixelToLatLng(long x, long y, int zoom)
221
      {
222
         PointLatLng ret = PointLatLng.Empty;
223
 
224
         var size = GetTileMatrixSizePixel(zoom);
225
 
226
         var oX = x << (20 - zoom);
227
         var oY = (size.Height - y) << (20 - zoom);
228
         {
229
            var l = PPToWGS(oX, oY);
230
            ret.Lat = Clip(l[0], MinLatitude, MaxLatitude);
231
            ret.Lng = Clip(l[1], MinLongitude, MaxLongitude);
232
         }
233
         return ret;
234
      }
235
 
236
      public override GSize GetTileMatrixSizeXY(int zoom)
237
      {
238
         return new GSize((long)Math.Pow(2, zoom), (long)Math.Pow(2, zoom));
239
      }
240
 
241
      public override GSize GetTileMatrixSizePixel(int zoom)
242
      {
243
         GSize s = GetTileMatrixSizeXY(zoom);
244
         return new GSize(s.Width << 8, s.Height << 8);
245
      }
246
 
247
      public override GSize GetTileMatrixMinXY(int zoom)
248
      {
249
         long wh = zoom > 3 ? (3 * (long)Math.Pow(2, zoom - 4)) : 1;
250
         return new GSize(wh, wh);
251
      }
252
 
253
      public override GSize GetTileMatrixMaxXY(int zoom)
254
      {
255
         long wh = (long)Math.Pow(2, zoom) - (long)Math.Pow(2, zoom - 2);
256
         return new GSize(wh, wh);
257
      }
258
   }
259
}