Subversion Repositories eFlore/Projets.eflore-projets

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
1001 delphine 1
<?php
2
/**
3
 * Author : Julien Moquet
4
 *
5
 * Inspired by Proj4js 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
/** datum object
11
 */
12
class proj4phpDatum {
13
 
14
    public $datum_type;
15
    public $datum_params;
16
 
17
    /**
18
     *
19
     * @param type $proj
20
     */
21
    public function __construct( $proj ) {
22
 
23
        $this->datum_type = Proj4php::$common->PJD_WGS84;   //default setting
24
 
25
        if( isset($proj->datumCode) && $proj->datumCode == 'none' ) {
26
            $this->datum_type = Proj4php::$common->PJD_NODATUM;
27
        }
28
 
29
        if( isset( $proj->datum_params ) ) {
30
 
31
            for( $i = 0; $i < sizeof( $proj->datum_params ); $i++ ) {
32
                $proj->datum_params[$i] = floatval( $proj->datum_params[$i] );
33
            }
34
 
35
            if( $proj->datum_params[0] != 0 || $proj->datum_params[1] != 0 || $proj->datum_params[2] != 0 ) {
36
                $this->datum_type = Proj4php::$common->PJD_3PARAM;
37
            }
38
 
39
            if( sizeof( $proj->datum_params ) > 3 ) {
40
                if( $proj->datum_params[3] != 0 || $proj->datum_params[4] != 0 ||
41
                    $proj->datum_params[5] != 0 || $proj->datum_params[6] != 0 ) {
42
 
43
                    $this->datum_type = Proj4php::$common->PJD_7PARAM;
44
                    $proj->datum_params[3] *= Proj4php::$common->SEC_TO_RAD;
45
                    $proj->datum_params[4] *= Proj4php::$common->SEC_TO_RAD;
46
                    $proj->datum_params[5] *= Proj4php::$common->SEC_TO_RAD;
47
                    $proj->datum_params[6] = ($proj->datum_params[6] / 1000000.0) + 1.0;
48
                }
49
            }
50
 
51
            $this->datum_params = $proj->datum_params;
52
        }
53
        if( isset( $proj ) ) {
54
            $this->a = $proj->a;    //datum object also uses these values
55
            $this->b = $proj->b;
56
            $this->es = $proj->es;
57
            $this->ep2 = $proj->ep2;
58
            #$this->datum_params = $proj->datum_params;
59
        }
60
    }
61
 
62
 
63
    /**
64
     *
65
     * @param type $dest
66
     * @return boolean Returns TRUE if the two datums match, otherwise FALSE.
67
     * @throws type
68
     */
69
    public function compare_datums( $dest ) {
70
 
71
        if( $this->datum_type != $dest->datum_type ) {
72
            return false; // false, datums are not equal
73
        } else if( $this->a != $dest->a || abs( $this->es - $dest->es ) > 0.000000000050 ) {
74
            // the tolerence for es is to ensure that GRS80 and WGS84
75
            // are considered identical
76
            return false;
77
        } else if( $this->datum_type == Proj4php::$common->PJD_3PARAM ) {
78
            return ($this->datum_params[0] == $dest->datum_params[0]
79
                      && $this->datum_params[1] == $dest->datum_params[1]
80
                      && $this->datum_params[2] == $dest->datum_params[2]);
81
        } else if( $this->datum_type == Proj4php::$common->PJD_7PARAM ) {
82
            return ($this->datum_params[0] == $dest->datum_params[0]
83
                      && $this->datum_params[1] == $dest->datum_params[1]
84
                      && $this->datum_params[2] == $dest->datum_params[2]
85
                      && $this->datum_params[3] == $dest->datum_params[3]
86
                      && $this->datum_params[4] == $dest->datum_params[4]
87
                      && $this->datum_params[5] == $dest->datum_params[5]
88
                      && $this->datum_params[6] == $dest->datum_params[6]);
89
        } else if( $this->datum_type == Proj4php::$common->PJD_GRIDSHIFT ||
90
            $dest->datum_type == Proj4php::$common->PJD_GRIDSHIFT ) {
91
            throw(new Exception( "ERROR: Grid shift transformations are not implemented." ));
92
            return false;
93
        }
94
 
95
        return true; // datums are equal
96
    }
97
 
98
    /*
99
     * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
100
     * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
101
     * according to the current ellipsoid parameters.
102
     *
103
     *    Latitude  : Geodetic latitude in radians                     (input)
104
     *    Longitude : Geodetic longitude in radians                    (input)
105
     *    Height    : Geodetic height, in meters                       (input)
106
     *    X         : Calculated Geocentric X coordinate, in meters    (output)
107
     *    Y         : Calculated Geocentric Y coordinate, in meters    (output)
108
     *    Z         : Calculated Geocentric Z coordinate, in meters    (output)
109
     *
110
     */
111
    public function geodetic_to_geocentric( $p ) {
112
 
113
        $Longitude = $p->x;
114
        $Latitude = $p->y;
115
        $Height = isset( $p->z ) ? $p->z : 0;   //Z value not always supplied
116
        $Error_Code = 0;  //  GEOCENT_NO_ERROR;
117
 
118
        /*
119
         * * Don't blow up if Latitude is just a little out of the value
120
         * * range as it may just be a rounding issue.  Also removed longitude
121
         * * test, it should be wrapped by cos() and sin().  NFW for PROJ.4, Sep/2001.
122
         */
123
        if( $Latitude < -Proj4php::$common->HALF_PI && $Latitude > -1.001 * Proj4php::$common->HALF_PI ) {
124
            $Latitude = -Proj4php::$common->HALF_PI;
125
        } else if( $Latitude > Proj4php::$common->HALF_PI && $Latitude < 1.001 * Proj4php::$common->HALF_PI ) {
126
            $Latitude = Proj4php::$common->HALF_PI;
127
        } else if( ($Latitude < -Proj4php::$common->HALF_PI) || ($Latitude > Proj4php::$common->HALF_PI) ) {
128
            /* Latitude out of range */
129
            Proj4php::reportError( 'geocent:lat out of range:' . $Latitude );
130
            return null;
131
        }
132
 
133
        if( $Longitude > Proj4php::$common->PI )
134
            $Longitude -= (2 * Proj4php::$common->PI);
135
 
136
        $Sin_Lat = sin( $Latitude ); /*  sin(Latitude)  */
137
        $Cos_Lat = cos( $Latitude ); /*  cos(Latitude)  */
138
        $Sin2_Lat = $Sin_Lat * $Sin_Lat; /*  Square of sin(Latitude)  */
139
        $Rn = $this->a / (sqrt( 1.0e0 - $this->es * $Sin2_Lat )); /*  Earth radius at location  */
140
        $p->x = ($Rn + $Height) * $Cos_Lat * cos( $Longitude );
141
        $p->y = ($Rn + $Height) * $Cos_Lat * sin( $Longitude );
142
        $p->z = (($Rn * (1 - $this->es)) + $Height) * $Sin_Lat;
143
 
144
        return $Error_Code;
145
    }
146
 
147
 
148
    /**
149
     *
150
     * @param object $p
151
     * @return type
152
     */
153
    public function geocentric_to_geodetic( $p ) {
154
 
155
        /* local defintions and variables */
156
        /* end-criterium of loop, accuracy of sin(Latitude) */
157
        $genau = 1.E-12;
158
        $genau2 = ($genau * $genau);
159
        $maxiter = 30;
160
        $X = $p->x;
161
        $Y = $p->y;
162
        $Z = $p->z ? $p->z : 0.0;   //Z value not always supplied
163
 
164
        /*
165
        $P;        // distance between semi-minor axis and location
166
        $RR;       // distance between center and location
167
        $CT;       // sin of geocentric latitude
168
        $ST;       // cos of geocentric latitude
169
        $RX;
170
        $RK;
171
        $RN;       // Earth radius at location
172
        $CPHI0;    // cos of start or old geodetic latitude in iterations
173
        $SPHI0;    // sin of start or old geodetic latitude in iterations
174
        $CPHI;     // cos of searched geodetic latitude
175
        $SPHI;     // sin of searched geodetic latitude
176
        $SDPHI;    // end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1))
177
        $At_Pole;     // indicates location is in polar region
178
        $iter;        // of continous iteration, max. 30 is always enough (s.a.)
179
        $Longitude;
180
        $Latitude;
181
        $Height;
182
        */
183
 
184
        $At_Pole = false;
185
        $P = sqrt( $X * $X + $Y * $Y );
186
        $RR = sqrt( $X * $X + $Y * $Y + $Z * $Z );
187
 
188
        /*      special cases for latitude and longitude */
189
        if( $P / $this->a < $genau ) {
190
 
191
            /*  special case, if P=0. (X=0., Y=0.) */
192
            $At_Pole = true;
193
            $Longitude = 0.0;
194
 
195
            /*  if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis
196
             *  of ellipsoid (=center of mass), Latitude becomes PI/2 */
197
            if( $RR / $this->a < $genau ) {
198
                $Latitude = Proj4php::$common->HALF_PI;
199
                $Height = -$this->b;
200
                return;
201
            }
202
        } else {
203
            /*  ellipsoidal (geodetic) longitude
204
             *  interval: -PI < Longitude <= +PI */
205
            $Longitude = atan2( $Y, $X );
206
        }
207
 
208
        /* --------------------------------------------------------------
209
         * Following iterative algorithm was developped by
210
         * "Institut für Erdmessung", University of Hannover, July 1988.
211
         * Internet: www.ife.uni-hannover.de
212
         * Iterative computation of CPHI,SPHI and Height.
213
         * Iteration of CPHI and SPHI to 10**-12 radian res$p->
214
         * 2*10**-7 arcsec.
215
         * --------------------------------------------------------------
216
         */
217
        $CT = $Z / $RR;
218
        $ST = $P / $RR;
219
        $RX = 1.0 / sqrt( 1.0 - $this->es * (2.0 - $this->es) * $ST * $ST );
220
        $CPHI0 = $ST * (1.0 - $this->es) * $RX;
221
        $SPHI0 = $CT * $RX;
222
        $iter = 0;
223
 
224
        /* loop to find sin(Latitude) res$p-> Latitude
225
         * until |sin(Latitude(iter)-Latitude(iter-1))| < genau */
226
        do {
227
            ++$iter;
228
            $RN = $this->a / sqrt( 1.0 - $this->es * $SPHI0 * $SPHI0 );
229
 
230
            /*  ellipsoidal (geodetic) height */
231
            $Height = $P * $CPHI0 + $Z * $SPHI0 - $RN * (1.0 - $this->es * $SPHI0 * $SPHI0);
232
 
233
            $RK = $this->es * $RN / ($RN + $Height);
234
            $RX = 1.0 / sqrt( 1.0 - $RK * (2.0 - $RK) * $ST * $ST );
235
            $CPHI = $ST * (1.0 - $RK) * $RX;
236
            $SPHI = $CT * $RX;
237
            $SDPHI = $SPHI * $CPHI0 - $CPHI * $SPHI0;
238
            $CPHI0 = $CPHI;
239
            $SPHI0 = $SPHI;
240
        } while( $SDPHI * $SDPHI > $genau2 && $iter < $maxiter );
241
 
242
        /*      ellipsoidal (geodetic) latitude */
243
        $Latitude = atan( $SPHI / abs( $CPHI ) );
244
 
245
        $p->x = $Longitude;
246
        $p->y = $Latitude;
247
        $p->z = $Height;
248
 
249
        return $p;
250
    }
251
 
252
    /**
253
     * Convert_Geocentric_To_Geodetic
254
     * The method used here is derived from 'An Improved Algorithm for
255
     * Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996
256
     *
257
     * @param object Point $p
258
     * @return object Point $p
259
     */
260
    public function geocentric_to_geodetic_noniter( $p ) {
261
 
262
        /*
263
        $Longitude;
264
        $Latitude;
265
        $Height;
266
 
267
        $W;        // distance from Z axis
268
        $W2;       // square of distance from Z axis
269
        $T0;       // initial estimate of vertical component
270
        $T1;       // corrected estimate of vertical component
271
        $S0;       // initial estimate of horizontal component
272
        $S1;       // corrected estimate of horizontal component
273
        $Sin_B0;   // sin(B0), B0 is estimate of Bowring aux variable
274
        $Sin3_B0;  // cube of sin(B0)
275
        $Cos_B0;   // cos(B0)
276
        $Sin_p1;   // sin(phi1), phi1 is estimated latitude
277
        $Cos_p1;   // cos(phi1)
278
        $Rn;       // Earth radius at location
279
        $Sum;      // numerator of cos(phi1)
280
        $At_Pole;  // indicates location is in polar region
281
        */
282
 
283
        $X = floatval( $p->x );  // cast from string to float
284
        $Y = floatval( $p->y );
285
        $Z = floatval( $p->z ? $p->z : 0 );
286
 
287
        $At_Pole = false;
288
        if( $X <> 0.0 ) {
289
            $Longitude = atan2( $Y, $X );
290
        } else {
291
            if( $Y > 0 ) {
292
                $Longitude = Proj4php::$common->HALF_PI;
293
            } else if( Y < 0 ) {
294
                $Longitude = -Proj4php::$common->HALF_PI;
295
            } else {
296
                $At_Pole = true;
297
                $Longitude = 0.0;
298
                if( $Z > 0.0 ) { /* north pole */
299
                    $Latitude = Proj4php::$common->HALF_PI;
300
                } else if( Z < 0.0 ) { /* south pole */
301
                    $Latitude = -Proj4php::$common->HALF_PI;
302
                } else { /* center of earth */
303
                    $Latitude = Proj4php::$common->HALF_PI;
304
                    $Height = -$this->b;
305
                    return;
306
                }
307
            }
308
        }
309
        $W2 = $X * $X + $Y * $Y;
310
        $W = sqrt( $W2 );
311
        $T0 = $Z * Proj4php::$common->AD_C;
312
        $S0 = sqrt( $T0 * $T0 + $W2 );
313
        $Sin_B0 = $T0 / $S0;
314
        $Cos_B0 = $W / $S0;
315
        $Sin3_B0 = $Sin_B0 * $Sin_B0 * $Sin_B0;
316
        $T1 = $Z + $this->b * $this->ep2 * $Sin3_B0;
317
        $Sum = $W - $this->a * $this->es * $Cos_B0 * $Cos_B0 * $Cos_B0;
318
        $S1 = sqrt( $T1 * $T1 + $Sum * $Sum );
319
        $Sin_p1 = $T1 / $S1;
320
        $Cos_p1 = $Sum / $S1;
321
        $Rn = $this->a / sqrt( 1.0 - $this->es * $Sin_p1 * $Sin_p1 );
322
        if( $Cos_p1 >= Proj4php::$common->COS_67P5 ) {
323
            $Height = $W / $Cos_p1 - $Rn;
324
        } else if( $Cos_p1 <= -Proj4php::$common->COS_67P5 ) {
325
            $Height = $W / -$Cos_p1 - $Rn;
326
        } else {
327
            $Height = $Z / $Sin_p1 + $Rn * ($this->es - 1.0);
328
        }
329
        if( $At_Pole == false ) {
330
            $Latitude = atan( $Sin_p1 / $Cos_p1 );
331
        }
332
 
333
        $p->x = $Longitude;
334
        $p->y = $Latitude;
335
        $p->z = $Height;
336
 
337
        return $p;
338
    }
339
 
340
    /************************************************************** */
341
    // pj_geocentic_to_wgs84( p )
342
    //  p = point to transform in geocentric coordinates (x,y,z)
343
    public function geocentric_to_wgs84( $p ) {
344
 
345
        if( $this->datum_type == Proj4php::$common->PJD_3PARAM ) {
346
            // if( x[io] == HUGE_VAL )
347
            //    continue;
348
            $p->x += $this->datum_params[0];
349
            $p->y += $this->datum_params[1];
350
            $p->z += $this->datum_params[2];
351
        } else if( $this->datum_type == Proj4php::$common->PJD_7PARAM ) {
352
            $Dx_BF = $this->datum_params[0];
353
            $Dy_BF = $this->datum_params[1];
354
            $Dz_BF = $this->datum_params[2];
355
            $Rx_BF = $this->datum_params[3];
356
            $Ry_BF = $this->datum_params[4];
357
            $Rz_BF = $this->datum_params[5];
358
            $M_BF = $this->datum_params[6];
359
            // if( x[io] == HUGE_VAL )
360
            //    continue;
361
            $p->x = $M_BF * ( $p->x - $Rz_BF * $p->y + $Ry_BF * $p->z) + $Dx_BF;
362
            $p->y = $M_BF * ( $Rz_BF * $p->x + $p->y - $Rx_BF * $p->z) + $Dy_BF;
363
            $p->z = $M_BF * (-$Ry_BF * $p->x + $Rx_BF * $p->y + $p->z) + $Dz_BF;
364
        }
365
    }
366
 
367
    /*************************************************************** */
368
 
369
    // pj_geocentic_from_wgs84()
370
    //  coordinate system definition,
371
    //  point to transform in geocentric coordinates (x,y,z)
372
    public function geocentric_from_wgs84( $p ) {
373
 
374
        if( $this->datum_type == Proj4php::$common->PJD_3PARAM ) {
375
            //if( x[io] == HUGE_VAL )
376
            //    continue;
377
            $p->x -= $this->datum_params[0];
378
            $p->y -= $this->datum_params[1];
379
            $p->z -= $this->datum_params[2];
380
        } else if( $this->datum_type == Proj4php::$common->PJD_7PARAM ) {
381
            $Dx_BF = $this->datum_params[0];
382
            $Dy_BF = $this->datum_params[1];
383
            $Dz_BF = $this->datum_params[2];
384
            $Rx_BF = $this->datum_params[3];
385
            $Ry_BF = $this->datum_params[4];
386
            $Rz_BF = $this->datum_params[5];
387
            $M_BF = $this->datum_params[6];
388
            $x_tmp = ($p->x - $Dx_BF) / $M_BF;
389
            $y_tmp = ($p->y - $Dy_BF) / $M_BF;
390
            $z_tmp = ($p->z - $Dz_BF) / $M_BF;
391
            //if( x[io] == HUGE_VAL )
392
            //    continue;
393
 
394
            $p->x = $x_tmp + $Rz_BF * $y_tmp - $Ry_BF * $z_tmp;
395
            $p->y = -$Rz_BF * $x_tmp + $y_tmp + $Rx_BF * $z_tmp;
396
            $p->z = $Ry_BF * $x_tmp - $Rx_BF * $y_tmp + $z_tmp;
397
        } //cs_geocentric_from_wgs84()
398
    }
399
 
400
}
401