Bienvenue sur PostGIS.fr

Bienvenue sur PostGIS.fr , le site de la communauté des utilisateurs francophones de PostGIS.

PostGIS ajoute le support d'objets géographique à la base de données PostgreSQL. En effet, PostGIS "spatialise" le serverur PostgreSQL, ce qui permet de l'utiliser comme une base de données SIG.

Maintenu à jour, en fonction de nos disponibilités et des diverses sorties des outils que nous testons, nous vous proposons l'ensemble de nos travaux publiés en langue française.

source: trunk/workshop-routing-foss4g/web/proj4js/lib/projCode/nzmg.js @ 76

Revision 76, 10.9 KB checked in by djay, 12 years ago (diff)

Ajout du répertoire web

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