Subversion Repositories eFlore/Projets.eflore-projets

Rev

Rev 1001 | Details | Compare with Previous | Last modification | View Log | RSS feed

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