| 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 |
// Initialize the Stereographic projection
|
|
|
11 |
class Proj4phpProjStere {
|
|
|
12 |
|
|
|
13 |
protected $TOL = 1.e-8;
|
|
|
14 |
protected $NITER = 8;
|
|
|
15 |
protected $CONV = 1.e-10;
|
|
|
16 |
protected $S_POLE = 0;
|
|
|
17 |
protected $N_POLE = 1;
|
|
|
18 |
protected $OBLIQ = 2;
|
|
|
19 |
protected $EQUIT = 3;
|
|
|
20 |
|
|
|
21 |
/**
|
|
|
22 |
*
|
|
|
23 |
* @param type $phit
|
|
|
24 |
* @param type $sinphi
|
|
|
25 |
* @param type $eccen
|
|
|
26 |
* @return type
|
|
|
27 |
*/
|
|
|
28 |
public function ssfn_( $phit, $sinphi, $eccen ) {
|
|
|
29 |
$sinphi *= $eccen;
|
|
|
30 |
return (tan( .5 * (Proj4php::$common->HALF_PI + $phit) ) * pow( (1. - $sinphi) / (1. + $sinphi), .5 * $eccen ));
|
|
|
31 |
}
|
|
|
32 |
|
|
|
33 |
/**
|
|
|
34 |
*
|
|
|
35 |
*/
|
|
|
36 |
public function init() {
|
|
|
37 |
$this->phits = $this->lat_ts ? $this->lat_ts : Proj4php::$common->HALF_PI;
|
|
|
38 |
$t = abs( $this->lat0 );
|
|
|
39 |
if( (abs( $t ) - Proj4php::$common->HALF_PI) < Proj4php::$common->EPSLN ) {
|
|
|
40 |
$this->mode = $this->lat0 < 0. ? $this->S_POLE : $this->N_POLE;
|
|
|
41 |
} else {
|
|
|
42 |
$this->mode = $t > Proj4php::$common->EPSLN ? $this->OBLIQ : $this->EQUIT;
|
|
|
43 |
}
|
|
|
44 |
$this->phits = abs( $this->phits );
|
|
|
45 |
if( $this->es ) {
|
|
|
46 |
#$X;
|
|
|
47 |
|
|
|
48 |
switch( $this->mode ) {
|
|
|
49 |
case $this->N_POLE:
|
|
|
50 |
case $this->S_POLE:
|
|
|
51 |
if( abs( $this->phits - Proj4php::$common->HALF_PI ) < Proj4php::$common->EPSLN ) {
|
|
|
52 |
$this->akm1 = 2. * $this->k0 / sqrt( pow( 1 + $this->e, 1 + $this->e ) * pow( 1 - $this->e, 1 - $this->e ) );
|
|
|
53 |
} else {
|
|
|
54 |
$t = sin( $this->phits );
|
|
|
55 |
$this->akm1 = cos( $this->phits ) / Proj4php::$common->tsfnz( $this->e, $this->phits, $t );
|
|
|
56 |
$t *= $this->e;
|
|
|
57 |
$this->akm1 /= sqrt( 1. - $t * $t );
|
|
|
58 |
}
|
|
|
59 |
break;
|
|
|
60 |
case $this->EQUIT:
|
|
|
61 |
$this->akm1 = 2. * $this->k0;
|
|
|
62 |
break;
|
|
|
63 |
case $this->OBLIQ:
|
|
|
64 |
$t = sin( $this->lat0 );
|
|
|
65 |
$X = 2. * atan( $this->ssfn_( $this->lat0, $t, $this->e ) ) - Proj4php::$common->HALF_PI;
|
|
|
66 |
$t *= $this->e;
|
|
|
67 |
$this->akm1 = 2. * $this->k0 * cos( $this->lat0 ) / sqrt( 1. - $t * $t );
|
|
|
68 |
$this->sinX1 = sin( $X );
|
|
|
69 |
$this->cosX1 = cos( $X );
|
|
|
70 |
break;
|
|
|
71 |
}
|
|
|
72 |
} else {
|
|
|
73 |
switch( $this->mode ) {
|
|
|
74 |
case $this->OBLIQ:
|
|
|
75 |
$this->sinph0 = sin( $this->lat0 );
|
|
|
76 |
$this->cosph0 = cos( $this->lat0 );
|
|
|
77 |
case $this->EQUIT:
|
|
|
78 |
$this->akm1 = 2. * $this->k0;
|
|
|
79 |
break;
|
|
|
80 |
case $this->S_POLE:
|
|
|
81 |
case $this->N_POLE:
|
|
|
82 |
$this->akm1 = abs( $this->phits - Proj4php::$common->HALF_PI ) >= Proj4php::$common->EPSLN ?
|
|
|
83 |
cos( $this->phits ) / tan( Proj4php::$common->FORTPI - .5 * $this->phits ) :
|
|
|
84 |
2. * $this->k0;
|
|
|
85 |
break;
|
|
|
86 |
}
|
|
|
87 |
}
|
|
|
88 |
}
|
|
|
89 |
|
|
|
90 |
/**
|
|
|
91 |
* Stereographic forward equations--mapping lat,long to x,y
|
|
|
92 |
*
|
|
|
93 |
* @param type $p
|
|
|
94 |
* @return type
|
|
|
95 |
*/
|
|
|
96 |
public function forward( $p ) {
|
|
|
97 |
|
|
|
98 |
$lon = $p->x;
|
|
|
99 |
$lon = Proj4php::$common->adjust_lon( $lon - $this->long0 );
|
|
|
100 |
$lat = $p->y;
|
|
|
101 |
#$x;
|
|
|
102 |
#$y;
|
|
|
103 |
|
|
|
104 |
if( $this->sphere ) {
|
|
|
105 |
/*
|
|
|
106 |
$sinphi;
|
|
|
107 |
$cosphi;
|
|
|
108 |
$coslam;
|
|
|
109 |
$sinlam;
|
|
|
110 |
*/
|
|
|
111 |
|
|
|
112 |
$sinphi = sin( $lat );
|
|
|
113 |
$cosphi = cos( $lat );
|
|
|
114 |
$coslam = cos( $lon );
|
|
|
115 |
$sinlam = sin( $lon );
|
|
|
116 |
switch( $this->mode ) {
|
|
|
117 |
case $this->EQUIT:
|
|
|
118 |
$y = 1. + $cosphi * $coslam;
|
|
|
119 |
if( y <= Proj4php::$common->EPSLN ) {
|
|
|
120 |
Proj4php::reportError("stere:forward:Equit");
|
|
|
121 |
}
|
|
|
122 |
$y = $this->akm1 / $y;
|
|
|
123 |
$x = $y * $cosphi * $sinlam;
|
|
|
124 |
$y *= $sinphi;
|
|
|
125 |
break;
|
|
|
126 |
case $this->OBLIQ:
|
|
|
127 |
$y = 1. + $this->sinph0 * $sinphi + $this->cosph0 * $cosphi * $coslam;
|
|
|
128 |
if( $y <= Proj4php::$common->EPSLN ) {
|
|
|
129 |
Proj4php::reportError("stere:forward:Obliq");
|
|
|
130 |
}
|
|
|
131 |
$y = $this->akm1 / $y;
|
|
|
132 |
$x = $y * $cosphi * $sinlam;
|
|
|
133 |
$y *= $this->cosph0 * $sinphi - $this->sinph0 * $cosphi * $coslam;
|
|
|
134 |
break;
|
|
|
135 |
case $this->N_POLE:
|
|
|
136 |
$coslam = -$coslam;
|
|
|
137 |
$lat = -$lat;
|
|
|
138 |
//Note no break here so it conitnues through S_POLE
|
|
|
139 |
case $this->S_POLE:
|
|
|
140 |
if( abs( $lat - Proj4php::$common->HALF_PI ) < $this->TOL ) {
|
|
|
141 |
Proj4php::reportError("stere:forward:S_POLE");
|
|
|
142 |
}
|
|
|
143 |
$y = $this->akm1 * tan( Proj4php::$common->FORTPI + .5 * $lat );
|
|
|
144 |
$x = $sinlam * $y;
|
|
|
145 |
$y *= $coslam;
|
|
|
146 |
break;
|
|
|
147 |
}
|
|
|
148 |
} else {
|
|
|
149 |
$coslam = cos( $lon );
|
|
|
150 |
$sinlam = sin( $lon );
|
|
|
151 |
$sinphi = sin( $lat );
|
|
|
152 |
if( $this->mode == $this->OBLIQ || $this->mode == $this->EQUIT ) {
|
|
|
153 |
$Xt = 2. * atan( $this->ssfn_( $lat, $sinphi, $this->e ) );
|
|
|
154 |
$sinX = sin( $Xt - Proj4php::$common->HALF_PI );
|
|
|
155 |
$cosX = cos( $Xt );
|
|
|
156 |
}
|
|
|
157 |
switch( $this->mode ) {
|
|
|
158 |
case $this->OBLIQ:
|
|
|
159 |
$A = $this->akm1 / ($this->cosX1 * (1. + $this->sinX1 * $sinX + $this->cosX1 * $cosX * $coslam));
|
|
|
160 |
$y = $A * ($this->cosX1 * $sinX - $this->sinX1 * $cosX * $coslam);
|
|
|
161 |
$x = $A * $cosX;
|
|
|
162 |
break;
|
|
|
163 |
case $this->EQUIT:
|
|
|
164 |
$A = 2. * $this->akm1 / (1. + $cosX * $coslam);
|
|
|
165 |
$y = $A * $sinX;
|
|
|
166 |
$x = $A * $cosX;
|
|
|
167 |
break;
|
|
|
168 |
case $this->S_POLE:
|
|
|
169 |
$lat = -$lat;
|
|
|
170 |
$coslam = - $coslam;
|
|
|
171 |
$sinphi = -$sinphi;
|
|
|
172 |
case $this->N_POLE:
|
|
|
173 |
$x = $this->akm1 * Proj4php::$common->tsfnz( $this->e, $lat, $sinphi );
|
|
|
174 |
$y = - $x * $coslam;
|
|
|
175 |
break;
|
|
|
176 |
}
|
|
|
177 |
$x = $x * $sinlam;
|
|
|
178 |
}
|
|
|
179 |
|
|
|
180 |
$p->x = $x * $this->a + $this->x0;
|
|
|
181 |
$p->y = $y * $this->a + $this->y0;
|
|
|
182 |
|
|
|
183 |
return $p;
|
|
|
184 |
}
|
|
|
185 |
|
|
|
186 |
|
|
|
187 |
/**
|
|
|
188 |
* Stereographic inverse equations--mapping x,y to lat/long
|
|
|
189 |
*
|
|
|
190 |
* @param type $p
|
|
|
191 |
* @return type
|
|
|
192 |
*/
|
|
|
193 |
public function inverse( $p ) {
|
|
|
194 |
$x = ($p->x - $this->x0) / $this->a; /* descale and de-offset */
|
|
|
195 |
$y = ($p->y - $this->y0) / $this->a;
|
|
|
196 |
/*
|
|
|
197 |
$lon;
|
|
|
198 |
$lat;
|
|
|
199 |
$cosphi;
|
|
|
200 |
$sinphi;
|
|
|
201 |
$rho;
|
|
|
202 |
$tp = 0.0;
|
|
|
203 |
$phi_l = 0.0;
|
|
|
204 |
$i;
|
|
|
205 |
*/
|
|
|
206 |
$halfe = 0.0;
|
|
|
207 |
$pi2 = 0.0;
|
|
|
208 |
|
|
|
209 |
if( $this->sphere ) {
|
|
|
210 |
/*
|
|
|
211 |
$c;
|
|
|
212 |
$rh;
|
|
|
213 |
$sinc;
|
|
|
214 |
$cosc;
|
|
|
215 |
*/
|
|
|
216 |
|
|
|
217 |
$rh = sqrt( $x * $x + $y * $y );
|
|
|
218 |
$c = 2. * atan( $rh / $this->akm1 );
|
|
|
219 |
$sinc = sin( $c );
|
|
|
220 |
$cosc = cos( $c );
|
|
|
221 |
$lon = 0.;
|
|
|
222 |
|
|
|
223 |
switch( $this->mode ) {
|
|
|
224 |
case $this->EQUIT:
|
|
|
225 |
if( abs( $rh ) <= Proj4php::$common->EPSLN ) {
|
|
|
226 |
$lat = 0.;
|
|
|
227 |
} else {
|
|
|
228 |
$lat = asin( $y * $sinc / $rh );
|
|
|
229 |
}
|
|
|
230 |
if( $cosc != 0. || $x != 0. )
|
|
|
231 |
$lon = atan2( $x * $sinc, $cosc * $rh );
|
|
|
232 |
break;
|
|
|
233 |
case $this->OBLIQ:
|
|
|
234 |
if( abs( $rh ) <= Proj4php::$common->EPSLN ) {
|
|
|
235 |
$lat = $this->phi0;
|
|
|
236 |
} else {
|
|
|
237 |
$lat = asin( $cosc * $this->sinph0 + $y * $sinc * $this->cosph0 / $rh );
|
|
|
238 |
}
|
|
|
239 |
$c = $cosc - $this->sinph0 * sin( $lat );
|
|
|
240 |
if( $c != 0. || $x != 0. ) {
|
|
|
241 |
$lon = atan2( $x * $sinc * $this->cosph0, $c * $rh );
|
|
|
242 |
}
|
|
|
243 |
break;
|
|
|
244 |
case $this->N_POLE:
|
|
|
245 |
$y = -$y;
|
|
|
246 |
case $this->S_POLE:
|
|
|
247 |
if( abs( $rh ) <= Proj4php::$common->EPSLN ) {
|
|
|
248 |
$lat = $this->phi0;
|
|
|
249 |
} else {
|
|
|
250 |
$lat = asin( $this->mode == $this->S_POLE ? -$cosc : $cosc );
|
|
|
251 |
}
|
|
|
252 |
$lon = ($x == 0. && $y == 0.) ? 0. : atan2( $x, $y );
|
|
|
253 |
break;
|
|
|
254 |
}
|
|
|
255 |
$p->x = Proj4php::$common->adjust_lon( $lon + $this->long0 );
|
|
|
256 |
$p->y = $lat;
|
|
|
257 |
} else {
|
|
|
258 |
$rho = sqrt( $x * $x + $y * $y );
|
|
|
259 |
switch( $this->mode ) {
|
|
|
260 |
case $this->OBLIQ:
|
|
|
261 |
case $this->EQUIT:
|
|
|
262 |
$tp = 2. * atan2( $rho * $this->cosX1, $this->akm1 );
|
|
|
263 |
$cosphi = cos( $tp );
|
|
|
264 |
$sinphi = sin( $tp );
|
|
|
265 |
if( $rho == 0.0 ) {
|
|
|
266 |
$phi_l = asin( $cosphi * $this->sinX1 );
|
|
|
267 |
} else {
|
|
|
268 |
$phi_l = asin( $cosphi * $this->sinX1 + ($y * $sinphi * $this->cosX1 / $rho) );
|
|
|
269 |
}
|
|
|
270 |
|
|
|
271 |
$tp = tan( .5 * (Proj4php::$common->HALF_PI + $phi_l) );
|
|
|
272 |
$x *= $sinphi;
|
|
|
273 |
$y = $rho * $this->cosX1 * $cosphi - $y * $this->sinX1 * $sinphi;
|
|
|
274 |
$pi2 = Proj4php::$common->HALF_PI;
|
|
|
275 |
$halfe = .5 * $this->e;
|
|
|
276 |
break;
|
|
|
277 |
case $this->N_POLE:
|
|
|
278 |
$y = -$y;
|
|
|
279 |
case $this->S_POLE:
|
|
|
280 |
$tp = - $rho / $this->akm1;
|
|
|
281 |
$phi_l = Proj4php::$common->HALF_PI - 2. * atan( $tp );
|
|
|
282 |
$pi2 = -Proj4php::$common->HALF_PI;
|
|
|
283 |
$halfe = -.5 * $this->e;
|
|
|
284 |
break;
|
|
|
285 |
}
|
|
|
286 |
for( $i = $this->NITER; $i--; $phi_l = $lat ) { //check this
|
|
|
287 |
$sinphi = $this->e * sin( $phi_l );
|
|
|
288 |
$lat = 2. * atan( $tp * pow( (1. + $sinphi) / (1. - $sinphi), $halfe ) ) - $pi2;
|
|
|
289 |
if( abs( phi_l - lat ) < $this->CONV ) {
|
|
|
290 |
if( $this->mode == $this->S_POLE )
|
|
|
291 |
$lat = -$lat;
|
|
|
292 |
$lon = ($x == 0. && $y == 0.) ? 0. : atan2( $x, $y );
|
|
|
293 |
$p->x = Proj4php::$common->adjust_lon( $lon + $this->long0 );
|
|
|
294 |
$p->y = $lat;
|
|
|
295 |
return $p;
|
|
|
296 |
}
|
|
|
297 |
}
|
|
|
298 |
}
|
|
|
299 |
}
|
|
|
300 |
|
|
|
301 |
}
|
|
|
302 |
|
|
|
303 |
Proj4php::$proj['stere'] = new Proj4phpProjStere();
|