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
  NAME                     ALBERS CONICAL EQUAL AREA
4
 
5
  PURPOSE:	Transforms input longitude and latitude to Easting and Northing
6
  for the Albers Conical Equal Area projection.  The longitude
7
  and latitude must be in radians.  The Easting and Northing
8
  values will be returned in meters.
9
 
10
  PROGRAMMER              DATE
11
  ----------              ----
12
  T. Mittan,       	Feb, 1992
13
 
14
  ALGORITHM REFERENCES
15
 
16
  1.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
17
  Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
18
  State Government Printing Office, Washington D.C., 1987.
19
 
20
  2.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
21
  U.S. Geological Survey Professional Paper 1453 , United State Government
22
  Printing Office, Washington D.C., 1989.
23
 *******************************************************************************/
24
 
25
/**
26
 * Author : Julien Moquet
27
 *
28
 * Inspired by Proj4php from Mike Adair madairATdmsolutions.ca
29
 *                      and Richard Greenwood rich@greenwoodma$p->com
30
 * License: LGPL as per: http://www.gnu.org/copyleft/lesser.html
31
 */
32
class Proj4phpProjAea {
33
 
34
    /**
35
     *
36
     * @return void
37
     */
38
    public function init() {
39
 
40
        if( abs( $this->lat1 + $this->lat2 ) < Proj4php::$common->EPSLN ) {
41
            Proj4php::reportError( "aeaInitEqualLatitudes" );
42
            return;
43
        }
44
        $this->temp = $this->b / $this->a;
45
        $this->es = 1.0 - pow( $this->temp, 2 );
46
        $this->e3 = sqrt( $this->es );
47
 
48
        $this->sin_po = sin( $this->lat1 );
49
        $this->cos_po = cos( $this->lat1 );
50
        $this->t1 = $this->sin_po;
51
        $this->con = $this->sin_po;
52
        $this->ms1 = Proj4php::$common->msfnz( $this->e3, $this->sin_po, $this->cos_po );
53
        $this->qs1 = Proj4php::$common->qsfnz( $this->e3, $this->sin_po, $this->cos_po );
54
 
55
        $this->sin_po = sin( $this->lat2 );
56
        $this->cos_po = cos( $this->lat2 );
57
        $this->t2 = $this->sin_po;
58
        $this->ms2 = Proj4php::$common->msfnz( $this->e3, $this->sin_po, $this->cos_po );
59
        $this->qs2 = Proj4php::$common->qsfnz( $this->e3, $this->sin_po, $this->cos_po );
60
 
61
        $this->sin_po = sin( $this->lat0 );
62
        $this->cos_po = cos( $this->lat0 );
63
        $this->t3 = $this->sin_po;
64
        $this->qs0 = Proj4php::$common->qsfnz( $this->e3, $this->sin_po, $this->cos_po );
65
 
66
        if( abs( $this->lat1 - $this->lat2 ) > Proj4php::$common->EPSLN ) {
67
            $this->ns0 = ($this->ms1 * $this->ms1 - $this->ms2 * $this->ms2) / ($this->qs2 - $this->qs1);
68
        } else {
69
            $this->ns0 = $this->con;
70
        }
71
 
72
        $this->c = $this->ms1 * $this->ms1 + $this->ns0 * $this->qs1;
73
        $this->rh = $this->a * sqrt( $this->c - $this->ns0 * $this->qs0 ) / $this->ns0;
74
    }
75
 
76
    /**
77
     * Albers Conical Equal Area forward equations--mapping lat,long to x,y
78
     *
79
     * @param Point $p
80
     * @return Point $p
81
     */
82
    public function forward( $p ) {
83
 
84
        $lon = $p->x;
85
        $lat = $p->y;
86
 
87
        $this->sin_phi = sin( $lat );
88
        $this->cos_phi = cos( $lat );
89
 
90
        $qs = Proj4php::$common->qsfnz( $this->e3, $this->sin_phi, $this->cos_phi );
91
        $rh1 = $this->a * sqrt( $this->c - $this->ns0 * $qs ) / $this->ns0;
92
        $theta = $this->ns0 * Proj4php::$common->adjust_lon( $lon - $this->long0 );
93
        $x = rh1 * sin( $theta ) + $this->x0;
94
        $y = $this->rh - $rh1 * cos( $theta ) + $this->y0;
95
 
96
        $p->x = $x;
97
        $p->y = $y;
98
 
99
        return $p;
100
    }
101
 
102
    /**
103
     *
104
     * @param Point $p
105
     * @return Point $p
106
     */
107
    public function inverse( $p ) {
108
 
109
        $p->x -= $this->x0;
110
        $p->y = $this->rh - $p->y + $this->y0;
111
 
112
        if( $this->ns0 >= 0 ) {
113
            $rh1 = sqrt( $p->x * $p->x + $p->y * $p->y );
114
            $con = 1.0;
115
        } else {
116
            $rh1 = -sqrt( $p->x * $p->x + $p->y * $p->y );
117
            $con = -1.0;
118
        }
119
 
120
        $theta = 0.0;
121
        if( $rh1 != 0.0 ) {
122
            $theta = atan2( $con * $p->x, $con * $p->y );
123
        }
124
 
125
        $con = $rh1 * $this->ns0 / $this->a;
126
        $qs = ($this->c - $con * $con) / $this->ns0;
127
 
128
        if( $this->e3 >= 1e-10 ) {
129
            $con = 1 - .5 * (1.0 - $this->es) * log( (1.0 - $this->e3) / (1.0 + $this->e3) ) / $this->e3;
130
            if( abs( abs( $con ) - abs( $qs ) ) > .0000000001 ) {
131
                $lat = $this->phi1z( $this->e3, $qs );
132
            } else {
133
                if( $qs >= 0 ) {
134
                    $lat = .5 * Proj4php::$Common->PI;
135
                } else {
136
                    $lat = -.5 * Proj4php::$Common->PI;
137
                }
138
            }
139
        } else {
140
            $lat = $this->phi1z( $this->e3, $qs );
141
        }
142
 
143
        $lon = Proj4php::$common->adjust_lon( $theta / $this->ns0 + $this->long0 );
144
 
145
        $p->x = $lon;
146
        $p->y = $lat;
147
 
148
        return $p;
149
    }
150
 
151
    /**
152
     * Function to compute phi1, the latitude for the inverse of the Albers Conical Equal-Area projection.
153
     *
154
     * @param type $eccent
155
     * @param type $qs
156
     * @return $phi or null on Convergence error
157
     */
158
    public function phi1z( $eccent, $qs ) {
159
 
160
        $phi = Proj4php::$common->asinz( .5 * $qs );
161
 
162
        if( $eccent < Proj4php::$common->EPSLN )
163
            return $phi;
164
 
165
        $eccnts = $eccent * $eccent;
166
        for( $i = 1; $i <= 25; ++$i ) {
167
            $sinphi = sin( $phi );
168
            $cosphi = cos( $phi );
169
            $con = $eccent * $sinphi;
170
            $com = 1.0 - $con * $con;
171
            $dphi = .5 * $com * $com / $cosphi * ($qs / (1.0 - $eccnts) - $sinphi / $com + .5 / $eccent * log( (1.0 - $con) / (1.0 + $con) ));
172
            $phi = $phi + $dphi;
173
            if( abs( $dphi ) <= 1e-7 )
174
                return $phi;
175
        }
176
 
177
        Proj4php::reportError( "aea:phi1z:Convergence error" );
178
 
179
        return null;
180
    }
181
 
182
}
183
 
184
Proj4php::$proj['aea'] = new Proj4phpProjAea();