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
/**
4
 * Author : Julien Moquet
5
 *
6
 * Inspired by Proj4php from Mike Adair madairATdmsolutions.ca
7
 *                      and Richard Greenwood rich@greenwoodma$p->com
8
 * License: LGPL as per: http://www.gnu.org/copyleft/lesser.html
9
 */
10
/* * *****************************************************************************
11
  NAME                  LAMBERT AZIMUTHAL EQUAL-AREA
12
 
13
  PURPOSE:	Transforms input longitude and latitude to Easting and
14
  Northing for the Lambert Azimuthal Equal-Area projection.  The
15
  longitude and latitude must be in radians.  The Easting
16
  and Northing values will be returned in meters.
17
 
18
  PROGRAMMER              DATE
19
  ----------              ----
20
  D. Steinwand, EROS      March, 1991
21
 
22
  This function was adapted from the Lambert Azimuthal Equal Area projection
23
  code (FORTRAN) in the General Cartographic Transformation Package software
24
  which is available from the U.S. Geological Survey National Mapping Division.
25
 
26
  ALGORITHM REFERENCES
27
 
28
  1.  "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder,
29
  The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355.
30
 
31
  2.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
32
  Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
33
  State Government Printing Office, Washington D.C., 1987.
34
 
35
  3.  "Software Documentation for GCTP General Cartographic Transformation
36
  Package", U.S. Geological Survey National Mapping Division, May 1982.
37
 * ***************************************************************************** */
38
 
39
class Proj4phpProjLaea {
40
 
41
    protected $S_POLE = 1;
42
    protected $N_POLE = 2;
43
    protected $EQUIT = 3;
44
    protected $OBLIQ = 4;
45
 
46
    protected $P00 = .33333333333333333333;
47
    protected $P01 = .17222222222222222222;
48
    protected $P02 = .10257936507936507936;
49
    protected $P10 = .06388888888888888888;
50
    protected $P11 = .06640211640211640211;
51
    protected $P20 = .01641501294219154443;
52
 
53
    /* Initialize the Lambert Azimuthal Equal Area projection
54
      ------------------------------------------------------ */
55
    public function init() {
56
        $t = abs( $this->lat0 );
57
        if( abs( $t - Proj4php::$common->HALF_PI ) < Proj4php::$common->EPSLN ) {
58
            $this->mode = $this->lat0 < 0. ? $this->S_POLE : $this->N_POLE;
59
        } else if( abs( $t ) < Proj4php::$common->EPSLN ) {
60
            $this->mode = $this->EQUIT;
61
        } else {
62
            $this->mode = $this->OBLIQ;
63
        }
64
        if( $this->es > 0 ) {
65
            #$sinphi;
66
 
67
            $this->qp = Proj4php::$common->qsfnz( $this->e, 1.0 );
68
            $this->mmf = .5 / (1. - $this->es);
69
            $this->apa = $this->authset( $this->es );
70
            switch( $this->mode ) {
71
                case $this->N_POLE:
72
                case $this->S_POLE:
73
                    $this->dd = 1.;
74
                    break;
75
                case $this->EQUIT:
76
                    $this->rq = sqrt( .5 * $this->qp );
77
                    $this->dd = 1. / $this->rq;
78
                    $this->xmf = 1.;
79
                    $this->ymf = .5 * $this->qp;
80
                    break;
81
                case $this->OBLIQ:
82
                    $this->rq = sqrt( .5 * $this->qp );
83
                    $sinphi = sin( $this->lat0 );
84
                    $this->sinb1 = Proj4php::$common->qsfnz( $this->e, $sinphi ) / $this->qp;
85
                    $this->cosb1 = sqrt( 1. - $this->sinb1 * $this->sinb1 );
86
                    $this->dd = cos( $this->lat0 ) / (sqrt( 1. - $this->es * $sinphi * $sinphi ) * $this->rq * $this->cosb1);
87
                    $this->ymf = ($this->xmf = $this->rq) / $this->dd;
88
                    $this->xmf *= $this->dd;
89
                    break;
90
            }
91
        } else {
92
            if( $this->mode == $this->OBLIQ ) {
93
                $this->sinph0 = sin( $this->lat0 );
94
                $this->cosph0 = cos( $this->lat0 );
95
            }
96
        }
97
    }
98
 
99
    /* Lambert Azimuthal Equal Area forward equations--mapping lat,long to x,y
100
      ----------------------------------------------------------------------- */
101
    public function forward( $p ) {
102
 
103
        /* Forward equations
104
          ----------------- */
105
        #$x;
106
        #$y;
107
        $lam = $p->x;
108
        $phi = $p->y;
109
        $lam = Proj4php::$common->adjust_lon( $lam - $this->long0 );
110
 
111
        if( $this->sphere ) {
112
            /*
113
            $coslam;
114
            $cosphi;
115
            $sinphi;
116
            */
117
 
118
            $sinphi = sin( $phi );
119
            $cosphi = cos( $phi );
120
            $coslam = cos( $lam );
121
            switch( $this->mode ) {
122
                case $this->OBLIQ:
123
                case $this->EQUIT:
124
                    $y = ($this->mode == $this->EQUIT) ? 1. + $cosphi * $coslam : 1. + $this->sinph0 * $sinphi + $this->cosph0 * $cosphi * $coslam;
125
                    if( y <= Proj4php::$common->EPSLN ) {
126
                        Proj4php::reportError( "laea:fwd:y less than eps" );
127
                        return null;
128
                    }
129
                    $y = sqrt( 2. / $y );
130
                    $x = $y * cosphi * sin( $lam );
131
                    $y *= ($this->mode == $this->EQUIT) ? $sinphi : $this->cosph0 * $sinphi - $this->sinph0 * $cosphi * $coslam;
132
                    break;
133
                case $this->N_POLE:
134
                    $coslam = -$coslam;
135
                case $this->S_POLE:
136
                    if( abs( $phi + $this->phi0 ) < Proj4php::$common->EPSLN ) {
137
                        Proj4php::reportError( "laea:fwd:phi < eps" );
138
                        return null;
139
                    }
140
                    $y = Proj4php::$common->FORTPI - $phi * .5;
141
                    $y = 2. * (($this->mode == $this->S_POLE) ? cos( $y ) : sin( $y ));
142
                    $x = $y * sin( $lam );
143
                    $y *= $coslam;
144
                    break;
145
            }
146
        } else {
147
            /*
148
            $coslam;
149
            $sinlam;
150
            $sinphi;
151
            $q;
152
            */
153
            $sinb = 0.0;
154
            $cosb = 0.0;
155
            $b = 0.0;
156
 
157
            $coslam = cos( $lam );
158
            $sinlam = sin( $lam );
159
            $sinphi = sin( $phi );
160
            $q = Proj4php::$common->qsfnz( $this->e, $sinphi );
161
            if( $this->mode == $this->OBLIQ || $this->mode == $this->EQUIT ) {
162
                $sinb = $q / $this->qp;
163
                $cosb = sqrt( 1. - $sinb * $sinb );
164
            }
165
            switch( $this->mode ) {
166
                case $this->OBLIQ:
167
                    $b = 1. + $this->sinb1 * $sinb + $this->cosb1 * $cosb * $coslam;
168
                    break;
169
                case $this->EQUIT:
170
                    $b = 1. + $cosb * $coslam;
171
                    break;
172
                case $this->N_POLE:
173
                    $b = Proj4php::$common->HALF_PI + $phi;
174
                    $q = $this->qp - $q;
175
                    break;
176
                case $this->S_POLE:
177
                    $b = $phi - Proj4php::$common->HALF_PI;
178
                    $q = $this->qp + $q;
179
                    break;
180
            }
181
            if( abs( $b ) < Proj4php::$common->EPSLN ) {
182
                Proj4php::reportError( "laea:fwd:b < eps" );
183
                return null;
184
            }
185
            switch( $this->mode ) {
186
                case $this->OBLIQ:
187
                case $this->EQUIT:
188
                    $b = sqrt( 2. / $b );
189
                    if( $this->mode == $this->OBLIQ ) {
190
                        $y = $this->ymf * $b * ($this->cosb1 * $sinb - $this->sinb1 * $cosb * $coslam);
191
                    } else {
192
                        $y = ($b = sqrt( 2. / (1. + $cosb * $coslam) )) * $sinb * $this->ymf;
193
                    }
194
                    $x = $this->xmf * $b * $cosb * $sinlam;
195
                    break;
196
                case $this->N_POLE:
197
                case $this->S_POLE:
198
                    if( q >= 0. ) {
199
                        $x = ($b = sqrt( $q )) * $sinlam;
200
                        $y = $coslam * (($this->mode == $this->S_POLE) ? $b : -$b);
201
                    } else {
202
                        $x = $y = 0.;
203
                    }
204
                    break;
205
            }
206
        }
207
 
208
        //v 1.0
209
        /*
210
          $sin_lat=sin(lat);
211
          $cos_lat=cos(lat);
212
 
213
          $sin_delta_lon=sin(delta_lon);
214
          $cos_delta_lon=cos(delta_lon);
215
 
216
          $g =$this->sin_lat_o * sin_lat +$this->cos_lat_o * cos_lat * cos_delta_lon;
217
          if (g == -1.0) {
218
          Proj4php::reportError("laea:fwd:Point projects to a circle of radius "+ 2.0 * R);
219
          return null;
220
          }
221
          $ksp = $this->a * sqrt(2.0 / (1.0 + g));
222
          $x = ksp * cos_lat * sin_delta_lon + $this->x0;
223
          $y = ksp * ($this->cos_lat_o * sin_lat - $this->sin_lat_o * cos_lat * cos_delta_lon) + $this->y0;
224
         */
225
        $p->x = $this->a * $x + $this->x0;
226
        $p->y = $this->a * $y + $this->y0;
227
 
228
        return $p;
229
    }
230
 
231
    /* Inverse equations
232
      ----------------- */
233
    public function inverse( $p ) {
234
        $p->x -= $this->x0;
235
        $p->y -= $this->y0;
236
        $x = $p->x / $this->a;
237
        $y = $p->y / $this->a;
238
 
239
        if( $this->sphere ) {
240
            $cosz = 0.0;
241
            #$rh;
242
            $sinz = 0.0;
243
 
244
            $rh = sqrt( $x * $x + $y * $y );
245
            $phi = $rh * .5;
246
            if( $phi > 1. ) {
247
                Proj4php::reportError( "laea:Inv:DataError" );
248
                return null;
249
            }
250
            $phi = 2. * asin( $phi );
251
            if( $this->mode == $this->OBLIQ || $this->mode == $this->EQUIT ) {
252
                $sinz = sin( $phi );
253
                $cosz = cos( $phi );
254
            }
255
            switch( $this->mode ) {
256
                case $this->EQUIT:
257
                    $phi = (abs( $rh ) <= Proj4php::$common->EPSLN) ? 0. : asin( $y * $sinz / $rh );
258
                    $x *= $sinz;
259
                    $y = $cosz * $rh;
260
                    break;
261
                case $this->OBLIQ:
262
                    $phi = (abs( $rh ) <= Proj4php::$common->EPSLN) ? $this->phi0 : asin( $cosz * $this->sinph0 + $y * $sinz * $this->cosph0 / $rh );
263
                    $x *= $sinz * $this->cosph0;
264
                    $y = ($cosz - sin( $phi ) * $this->sinph0) * $rh;
265
                    break;
266
                case $this->N_POLE:
267
                    $y = -$y;
268
                    $phi = Proj4php::$common->HALF_PI - $phi;
269
                    break;
270
                case $this->S_POLE:
271
                    $phi -= Proj4php::$common->HALF_PI;
272
                    break;
273
            }
274
            $lam = ($y == 0. && ($this->mode == $this->EQUIT || $this->mode == $this->OBLIQ)) ? 0. : atan2( $x, $y );
275
        } else {
276
            /*
277
            $cCe;
278
            $sCe;
279
            $q;
280
            $rho;
281
            */
282
            $ab = 0.0;
283
 
284
            switch( $this->mode ) {
285
                case $this->EQUIT:
286
                case $this->OBLIQ:
287
                    $x /= $this->dd;
288
                    $y *= $this->dd;
289
                    $rho = sqrt( $x * $x + $y * $y );
290
                    if( $rho < Proj4php::$common->EPSLN ) {
291
                        $p->x = 0.;
292
                        $p->y = $this->phi0;
293
                        return $p;
294
                    }
295
                    $sCe = 2. * asin( .5 * $rho / $this->rq );
296
                    $cCe = cos( $sCe );
297
                    $x *= ($sCe = sin( $sCe ));
298
                    if( $this->mode == $this->OBLIQ ) {
299
                        $ab = $cCe * $this->sinb1 + $y * $sCe * $this->cosb1 / $rho;
300
                        $q = $this->qp * $ab;
301
                        $y = $rho * $this->cosb1 * $cCe - $y * $this->sinb1 * $sCe;
302
                    } else {
303
                        $ab = $y * $sCe / $rho;
304
                        $q = $this->qp * $ab;
305
                        $y = $rho * $cCe;
306
                    }
307
                    break;
308
                case $this->N_POLE:
309
                    $y = -$y;
310
                case $this->S_POLE:
311
                    $q = ($x * $x + $y * $y);
312
                    if( !$q ) {
313
                        $p->x = 0.;
314
                        $p->y = $this->phi0;
315
                        return $p;
316
                    }
317
                    /*
318
                      q = $this->qp - q;
319
                     */
320
                    $ab = 1. - $q / $this->qp;
321
                    if( $this->mode == $this->S_POLE ) {
322
                        $ab = - $ab;
323
                    }
324
                    break;
325
            }
326
            $lam = atan2( $x, $y );
327
            $phi = $this->authlat( asin( $ab ), $this->apa );
328
        }
329
 
330
        /*
331
          $Rh = sqrt($p->x *$p->x +$p->y * $p->y);
332
          $temp = Rh / (2.0 * $this->a);
333
 
334
          if (temp > 1) {
335
          Proj4php::reportError("laea:Inv:DataError");
336
          return null;
337
          }
338
 
339
          $z = 2.0 * Proj4php::$common.asinz(temp);
340
          $sin_z=sin(z);
341
          $cos_z=cos(z);
342
 
343
          $lon =$this->long0;
344
          if (abs(Rh) > Proj4php::$common->EPSLN) {
345
          $lat = Proj4php::$common.asinz($this->sin_lat_o * cos_z +$this-> cos_lat_o * sin_z *$p->y / Rh);
346
          $temp =abs($this->lat0) - Proj4php::$common->HALF_PI;
347
          if (abs(temp) > Proj4php::$common->EPSLN) {
348
          temp = cos_z -$this->sin_lat_o * sin(lat);
349
          if(temp!=0.0) lon=Proj4php::$common->adjust_lon($this->long0+atan2($p->x*sin_z*$this->cos_lat_o,temp*Rh));
350
          } else if ($this->lat0 < 0.0) {
351
          lon = Proj4php::$common->adjust_lon($this->long0 - atan2(-$p->x,$p->y));
352
          } else {
353
          lon = Proj4php::$common->adjust_lon($this->long0 + atan2($p->x, -$p->y));
354
          }
355
          } else {
356
          lat = $this->lat0;
357
          }
358
         */
359
        //return(OK);
360
        $p->x = Proj4php::$common->adjust_lon( $this->long0 + $lam );
361
        $p->y = $phi;
362
        return $p;
363
    }
364
 
365
    /**
366
     * determine latitude from authalic latitude
367
     *
368
     * @param type $es
369
     * @return type
370
     */
371
    public function authset( $es ) {
372
        #$t;
373
        $APA = array( );
374
        $APA[0] = $es * $this->P00;
375
        $t = $es * $es;
376
        $APA[0] += $t * $this->P01;
377
        $APA[1] = $t * $this->P10;
378
        $t *= $es;
379
        $APA[0] += $t * $this->P02;
380
        $APA[1] += $t * $this->P11;
381
        $APA[2] = $t * $this->P20;
382
        return $APA;
383
    }
384
 
385
    /**
386
     *
387
     * @param type $beta
388
     * @param type $APA
389
     * @return type
390
     */
391
    public function authlat( $beta, $APA ) {
392
        $t = $beta + $beta;
393
        return($beta + $APA[0] * sin( $t ) + $APA[1] * sin( $t + $t ) + $APA[2] * sin( $t + $t + $t ));
394
    }
395
 
396
}
397
 
398
Proj4php::$proj['laea'] = new Proj4phpProjLaea();