Understand shapefiles in different projections (notably OSGB)
[potlatch2.git] / com / gradoservice / proj4as / proj / ProjNzmg.as
1 /*******************************************************************************\r
2 NAME                            NEW ZEALAND MAP GRID\r
3 \r
4 PURPOSE:        Transforms input longitude and latitude to Easting and\r
5                 Northing for the New Zealand Map Grid projection.  The\r
6                 longitude and latitude must be in radians.  The Easting\r
7                 and Northing values will be returned in meters.\r
8 \r
9 \r
10 ALGORITHM REFERENCES\r
11 \r
12 1.  Department of Land and Survey Technical Circular 1973/32\r
13       http://www.linz.govt.nz/docs/miscellaneous/nz-map-definition.pdf\r
14 \r
15 2.  OSG Technical Report 4.1\r
16       http://www.linz.govt.nz/docs/miscellaneous/nzmg.pdf\r
17 \r
18 \r
19 IMPLEMENTATION NOTES\r
20 \r
21 The two references use different symbols for the calculated values. This\r
22 implementation uses the variable names similar to the symbols in reference [1].\r
23 \r
24 The alogrithm uses different units for delta latitude and delta longitude.\r
25 The delta latitude is assumed to be in units of seconds of arc x 10^-5.\r
26 The delta longitude is the usual radians. Look out for these conversions.\r
27 \r
28 The algorithm is described using complex arithmetic. There were three\r
29 options:\r
30    * find and use a Javascript library for complex arithmetic\r
31    * write my own complex library\r
32    * expand the complex arithmetic by hand to simple arithmetic\r
33 \r
34 This implementation has expanded the complex multiplication operations\r
35 into parallel simple arithmetic operations for the real and imaginary parts.\r
36 The imaginary part is way over to the right of the display; this probably\r
37 violates every coding standard in the world, but, to me, it makes it much\r
38 more obvious what is going on.\r
39 \r
40 The following complex operations are used:\r
41    - addition\r
42    - multiplication\r
43    - division\r
44    - complex number raised to integer power\r
45    - summation\r
46 \r
47 A summary of complex arithmetic operations:\r
48    (from http://en.wikipedia.org/wiki/Complex_arithmetic)\r
49    addition:       (a + bi) + (c + di) = (a + c) + (b + d)i\r
50    subtraction:    (a + bi) - (c + di) = (a - c) + (b - d)i\r
51    multiplication: (a + bi) x (c + di) = (ac - bd) + (bc + ad)i\r
52    division:       (a + bi) / (c + di) = [(ac + bd)/(cc + dd)] + [(bc - ad)/(cc + dd)]i\r
53 \r
54 The algorithm needs to calculate summations of simple and complex numbers. This is\r
55 implemented using a for-loop, pre-loading the summed value to zero.\r
56 \r
57 The algorithm needs to calculate theta^2, theta^3, etc while doing a summation.\r
58 There are three possible implementations:\r
59    - use Math.pow in the summation loop - except for complex numbers\r
60    - precalculate the values before running the loop\r
61    - calculate theta^n = theta^(n-1) * theta during the loop\r
62 This implementation uses the third option for both real and complex arithmetic.\r
63 \r
64 For example\r
65    psi_n = 1;\r
66    sum = 0;\r
67    for (n = 1; n <=6; n++) {\r
68       psi_n1 = psi_n * psi;       // calculate psi^(n+1)\r
69       psi_n = psi_n1;\r
70       sum = sum + A[n] * psi_n;\r
71    }\r
72 \r
73 \r
74 TEST VECTORS\r
75 \r
76 NZMG E, N:         2487100.638      6751049.719     metres\r
77 NZGD49 long, lat:      172.739194       -34.444066  degrees\r
78 \r
79 NZMG E, N:         2486533.395      6077263.661     metres\r
80 NZGD49 long, lat:      172.723106       -40.512409  degrees\r
81 \r
82 NZMG E, N:         2216746.425      5388508.765     metres\r
83 NZGD49 long, lat:      169.172062       -46.651295  degrees\r
84 \r
85 Note that these test vectors convert from NZMG metres to lat/long referenced\r
86 to NZGD49, not the more usual WGS84. The difference is about 70m N/S and about\r
87 10m E/W.\r
88 \r
89 These test vectors are provided in reference [1]. Many more test\r
90 vectors are available in\r
91    http://www.linz.govt.nz/docs/topography/topographicdata/placenamesdatabase/nznamesmar08.zip\r
92 which is a catalog of names on the 260-series maps.\r
93 \r
94 \r
95 EPSG CODES\r
96 \r
97 NZMG     EPSG:27200\r
98 NZGD49   EPSG:4272\r
99 \r
100 http://spatialreference.org/ defines these as\r
101   Proj4js.defs["EPSG:4272"] = "+proj=longlat +ellps=intl +datum=nzgd49 +no_defs ";\r
102   Proj4js.defs["EPSG:27200"] = "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +ellps=intl +datum=nzgd49 +units=m +no_defs ";\r
103 \r
104 \r
105 LICENSE\r
106   Copyright: Stephen Irons 2008\r
107   Released under terms of the LGPL as per: http://www.gnu.org/copyleft/lesser.html\r
108 \r
109 *******************************************************************************/\r
110 \r
111 package com.gradoservice.proj4as.proj\r
112 {\r
113         import com.gradoservice.proj4as.ProjPoint;\r
114         import com.gradoservice.proj4as.ProjConstants;\r
115         import com.gradoservice.proj4as.Datum;\r
116                 \r
117         public class ProjNzmg extends AbstractProjProjection\r
118         {\r
119                 private var A:Array = [];\r
120                 private var B_re:Array = [];\r
121                 private var B_im:Array = [];\r
122                 private var C_re:Array = [];\r
123                 private var C_im:Array = [];\r
124                 private var D:Array = [];\r
125                 private var iterations:int = 1;                                         \r
126                 \r
127                 public function ProjNzmg(data:ProjParams)\r
128                 {\r
129                         super(data);\r
130                 }\r
131                 \r
132                   override public function init():void\r
133                   {\r
134                     this.A[1]  = +0.6399175073;\r
135                     this.A[2]  = -0.1358797613;\r
136                     this.A[3]  = +0.063294409;\r
137                     this.A[4]  = -0.02526853;\r
138                     this.A[5]  = +0.0117879;\r
139                     this.A[6]  = -0.0055161;\r
140                     this.A[7]  = +0.0026906;\r
141                     this.A[8]  = -0.001333;\r
142                     this.A[9]  = +0.00067;\r
143                     this.A[10] = -0.00034;\r
144                 \r
145                     this.B_re[1] = +0.7557853228;   this.B_im[1] =  0.0;\r
146                     this.B_re[2] = +0.249204646;    this.B_im[2] = +0.003371507;\r
147                     this.B_re[3] = -0.001541739;    this.B_im[3] = +0.041058560;\r
148                     this.B_re[4] = -0.10162907;     this.B_im[4] = +0.01727609;\r
149                     this.B_re[5] = -0.26623489;     this.B_im[5] = -0.36249218;\r
150                     this.B_re[6] = -0.6870983;      this.B_im[6] = -1.1651967;\r
151                 \r
152                     this.C_re[1] = +1.3231270439;   this.C_im[1] =  0.0;\r
153                     this.C_re[2] = -0.577245789;    this.C_im[2] = -0.007809598;\r
154                     this.C_re[3] = +0.508307513;    this.C_im[3] = -0.112208952;\r
155                     this.C_re[4] = -0.15094762;     this.C_im[4] = +0.18200602;\r
156                     this.C_re[5] = +1.01418179;     this.C_im[5] = +1.64497696;\r
157                     this.C_re[6] = +1.9660549;      this.C_im[6] = +2.5127645;\r
158                 \r
159                     this.D[1] = +1.5627014243;\r
160                     this.D[2] = +0.5185406398;\r
161                     this.D[3] = -0.03333098;\r
162                     this.D[4] = -0.1052906;\r
163                     this.D[5] = -0.0368594;\r
164                     this.D[6] = +0.007317;\r
165                     this.D[7] = +0.01220;\r
166                     this.D[8] = +0.00394;\r
167                     this.D[9] = -0.0013;\r
168                   }\r
169                 \r
170                   /**\r
171                     New Zealand Map Grid Forward  - long/lat to x/y\r
172                     long/lat in radians\r
173                   */\r
174                   override public function forward(p:ProjPoint):ProjPoint\r
175                   {\r
176                     var lon:Number = p.x;\r
177                     var lat:Number = p.y;\r
178                 \r
179                     var delta_lat:Number = lat - this.lat0;\r
180                     var delta_lon:Number = lon - this.long0;\r
181                 \r
182                     // 1. Calculate d_phi and d_psi    ...                          // and d_lambda\r
183                     // For this algorithm, delta_latitude is in seconds of arc x 10-5, so we need to scale to those units. Longitude is radians.\r
184                     var d_phi:Number = delta_lat / ProjConstants.SEC_TO_RAD * 1E-5;       var d_lambda:Number = delta_lon;\r
185                     var d_phi_n:Number = 1;  // d_phi^0\r
186                 \r
187                     var d_psi:Number = 0;\r
188                     for (var n:int = 1; n <= 10; n++) {\r
189                       d_phi_n = d_phi_n * d_phi;\r
190                       d_psi = d_psi + this.A[n] * d_phi_n;\r
191                     }\r
192                 \r
193                     // 2. Calculate theta\r
194                     var th_re:Number = d_psi;                                              var th_im:Number = d_lambda;\r
195                 \r
196                     // 3. Calculate z\r
197                     var th_n_re:Number = 1;                                                var th_n_im:Number = 0;  // theta^0\r
198                     var th_n_re1:Number;                                                   var th_n_im1:Number;\r
199                 \r
200                     var z_re:Number = 0;                                                   var z_im:Number = 0;\r
201                     for (n = 1; n <= 6; n++) {\r
202                       th_n_re1 = th_n_re*th_re - th_n_im*th_im;                     th_n_im1 = th_n_im*th_re + th_n_re*th_im;\r
203                       th_n_re = th_n_re1;                                           th_n_im = th_n_im1;\r
204                       z_re = z_re + this.B_re[n]*th_n_re - this.B_im[n]*th_n_im;    z_im = z_im + this.B_im[n]*th_n_re + this.B_re[n]*th_n_im;\r
205                     }\r
206                 \r
207                     // 4. Calculate easting and northing\r
208                     var x:Number = (z_im * this.a) + this.x0;\r
209                     var y:Number = (z_re * this.a) + this.y0;\r
210                 \r
211                     p.x = x; p.y = y;\r
212                 \r
213                     return p;\r
214                   }\r
215                 \r
216                 \r
217                   /**\r
218                     New Zealand Map Grid Inverse  -  x/y to long/lat\r
219                   */\r
220                  override public function inverse(p:ProjPoint):ProjPoint\r
221                  {\r
222                     var x:Number = p.x;\r
223                     var y:Number = p.y;\r
224                 \r
225                     var delta_x:Number = x - this.x0;\r
226                     var delta_y:Number = y - this.y0;\r
227                 \r
228                     // 1. Calculate z\r
229                     var z_re:Number = delta_y / this.a;                                              var z_im:Number = delta_x / this.a;\r
230                 \r
231                     // 2a. Calculate theta - first approximation gives km accuracy\r
232                     var z_n_re:Number = 1;                                                           var z_n_im:Number = 0;  // z^0\r
233                     var z_n_re1:Number;                                                              var z_n_im1:Number;\r
234                 \r
235                     var th_re:Number = 0;                                                            var th_im:Number = 0;\r
236                     for (var n:int = 1; n <= 6; n++) {\r
237                       z_n_re1 = z_n_re*z_re - z_n_im*z_im;                                    z_n_im1 = z_n_im*z_re + z_n_re*z_im;\r
238                       z_n_re = z_n_re1;                                                       z_n_im = z_n_im1;\r
239                       th_re = th_re + this.C_re[n]*z_n_re - this.C_im[n]*z_n_im;              th_im = th_im + this.C_im[n]*z_n_re + this.C_re[n]*z_n_im;\r
240                     }\r
241                 \r
242                     // 2b. Iterate to refine the accuracy of the calculation\r
243                     //        0 iterations gives km accuracy\r
244                     //        1 iteration gives m accuracy -- good enough for most mapping applications\r
245                     //        2 iterations bives mm accuracy\r
246                     for (var i:int = 0; i < this.iterations; i++) {\r
247                        var th_n_re:Number = th_re;                                                      var th_n_im:Number = th_im;\r
248                        var th_n_re1:Number;                                                             var th_n_im1:Number;\r
249                 \r
250                        var num_re:Number = z_re;                                                        var num_im:Number = z_im;\r
251                        for (n = 2; n <= 6; n++) {\r
252                          th_n_re1 = th_n_re*th_re - th_n_im*th_im;                               th_n_im1 = th_n_im*th_re + th_n_re*th_im;\r
253                          th_n_re = th_n_re1;                                                     th_n_im = th_n_im1;\r
254                          num_re = num_re + (n-1)*(this.B_re[n]*th_n_re - this.B_im[n]*th_n_im);  num_im = num_im + (n-1)*(this.B_im[n]*th_n_re + this.B_re[n]*th_n_im);\r
255                        }\r
256                 \r
257                        th_n_re = 1;                                                              th_n_im = 0;\r
258                        var den_re:Number = this.B_re[1];                                                var den_im:Number = this.B_im[1];\r
259                        for (n = 2; n <= 6; n++) {\r
260                          th_n_re1 = th_n_re*th_re - th_n_im*th_im;                               th_n_im1 = th_n_im*th_re + th_n_re*th_im;\r
261                          th_n_re = th_n_re1;                                                     th_n_im = th_n_im1;\r
262                          den_re = den_re + n * (this.B_re[n]*th_n_re - this.B_im[n]*th_n_im);    den_im = den_im + n * (this.B_im[n]*th_n_re + this.B_re[n]*th_n_im);\r
263                        }\r
264                 \r
265                        // Complex division\r
266                        var den2:Number = den_re*den_re + den_im*den_im;\r
267                        th_re = (num_re*den_re + num_im*den_im) / den2;                           th_im = (num_im*den_re - num_re*den_im) / den2;\r
268                     }\r
269                 \r
270                     // 3. Calculate d_phi              ...                                    // and d_lambda\r
271                     var d_psi:Number = th_re;                                                        var d_lambda:Number = th_im;\r
272                     var d_psi_n:Number = 1;  // d_psi^0\r
273                 \r
274                     var d_phi:Number = 0;\r
275                     for (n = 1; n <= 9; n++) {\r
276                        d_psi_n = d_psi_n * d_psi;\r
277                        d_phi = d_phi + this.D[n] * d_psi_n;\r
278                     }\r
279                 \r
280                     // 4. Calculate latitude and longitude\r
281                     // d_phi is calcuated in second of arc * 10^-5, so we need to scale back to radians. d_lambda is in radians.\r
282                     var lat:Number = this.lat0 + (d_phi * ProjConstants.SEC_TO_RAD * 1E5);\r
283                     var lon:Number = this.long0 +  d_lambda;\r
284                 \r
285                     p.x = lon;\r
286                     p.y = lat;\r
287                 \r
288                     return p;\r
289                   }             \r
290                 \r
291         }\r
292 }