| 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 SWISS OBLIQUE MERCATOR
|
|
|
11 |
|
|
|
12 |
PURPOSE: Swiss projection.
|
|
|
13 |
WARNING: X and Y are inverted (weird) in the swiss coordinate system. Not
|
|
|
14 |
here, since we want X to be horizontal and Y vertical.
|
|
|
15 |
|
|
|
16 |
ALGORITHM REFERENCES
|
|
|
17 |
1. "Formules et constantes pour le Calcul pour la
|
|
|
18 |
projection cylindrique conforme à axe oblique et pour la transformation entre
|
|
|
19 |
des systèmes de référence".
|
|
|
20 |
http://www.swisstopo.admin.ch/internet/swisstopo/fr/home/topics/survey/sys/refsys/switzerland.parsysrelated1.31216.downloadList.77004.DownloadFile.tmp/swissprojectionfr.pdf
|
|
|
21 |
|
|
|
22 |
*******************************************************************************/
|
|
|
23 |
|
|
|
24 |
class Proj4phpProjSomerc {
|
|
|
25 |
|
|
|
26 |
/**
|
|
|
27 |
*
|
|
|
28 |
*/
|
|
|
29 |
public function init() {
|
|
|
30 |
$phy0 = $this->lat0;
|
|
|
31 |
$this->lambda0 = $this->long0;
|
|
|
32 |
$sinPhy0 = sin( $phy0 );
|
|
|
33 |
$semiMajorAxis = $this->a;
|
|
|
34 |
$invF = $this->rf;
|
|
|
35 |
$flattening = 1 / $invF;
|
|
|
36 |
$e2 = 2 * $flattening - pow( $flattening, 2 );
|
|
|
37 |
$e = $this->e = sqrt( $e2 );
|
|
|
38 |
$this->R = $this->k0 * $semiMajorAxis * sqrt( 1 - $e2 ) / (1 - $e2 * pow( $sinPhy0, 2.0 ));
|
|
|
39 |
$this->alpha = sqrt( 1 + $e2 / (1 - $e2) * pow( cos( $phy0 ), 4.0 ) );
|
|
|
40 |
$this->b0 = asin( $sinPhy0 / $this->alpha );
|
|
|
41 |
$this->K = log( tan( $PI / 4.0 + $this->b0 / 2.0 ) )
|
|
|
42 |
- $this->alpha
|
|
|
43 |
* log( tan( $PI / 4.0 + $phy0 / 2.0 ) )
|
|
|
44 |
+ $this->alpha
|
|
|
45 |
* $e / 2
|
|
|
46 |
* log( (1 + $e * $sinPhy0)
|
|
|
47 |
/ (1 - $e * $sinPhy0) );
|
|
|
48 |
}
|
|
|
49 |
|
|
|
50 |
/**
|
|
|
51 |
*
|
|
|
52 |
* @param type $p
|
|
|
53 |
* @return type
|
|
|
54 |
*/
|
|
|
55 |
public function forward( $p ) {
|
|
|
56 |
$Sa1 = log( tan( $PI / 4.0 - $p->y / 2.0 ) );
|
|
|
57 |
$Sa2 = $this->e / 2.0
|
|
|
58 |
* log( (1 + $this->e * sin( $p->y ))
|
|
|
59 |
/ (1 - $this->e * sin( $p->y )) );
|
|
|
60 |
$S = -$this->alpha * ($Sa1 + $Sa2) + $this->K;
|
|
|
61 |
|
|
|
62 |
// spheric latitude
|
|
|
63 |
$b = 2.0 * (atan( exp( $S ) ) - proj4phpCommon::PI / 4.0);
|
|
|
64 |
|
|
|
65 |
// spheric longitude
|
|
|
66 |
$I = $this->alpha * ($p->x - $this->lambda0);
|
|
|
67 |
|
|
|
68 |
// psoeudo equatorial rotation
|
|
|
69 |
$rotI = atan( sin( $I )
|
|
|
70 |
/ (sin( $this->b0 ) * tan( $b ) +
|
|
|
71 |
cos( $this->b0 ) * cos( $I )) );
|
|
|
72 |
|
|
|
73 |
$rotB = asin( cos( $this->b0 ) * sin( $b ) -
|
|
|
74 |
sin( $this->b0 ) * cos( $b ) * cos( $I ) );
|
|
|
75 |
|
|
|
76 |
$p->y = $this->R / 2.0
|
|
|
77 |
* log( (1 + sin( $rotB )) / (1 - sin( $rotB )) )
|
|
|
78 |
+ $this->y0;
|
|
|
79 |
|
|
|
80 |
$p->x = $this->R * $rotI + $this->x0;
|
|
|
81 |
|
|
|
82 |
return $p;
|
|
|
83 |
}
|
|
|
84 |
|
|
|
85 |
/**
|
|
|
86 |
*
|
|
|
87 |
* @param type $p
|
|
|
88 |
* @return type
|
|
|
89 |
*/
|
|
|
90 |
public function inverse( $p ) {
|
|
|
91 |
|
|
|
92 |
$Y = $p->x - $this->x0;
|
|
|
93 |
$X = $p->y - $this->y0;
|
|
|
94 |
|
|
|
95 |
$rotI = $Y / $this->R;
|
|
|
96 |
$rotB = 2 * (atan( exp( $X / $this->R ) ) - $PI / 4.0);
|
|
|
97 |
|
|
|
98 |
$b = asin( cos( $this->b0 ) * sin( $rotB )
|
|
|
99 |
+ sin( $this->b0 ) * cos( $rotB ) * cos( $rotI ) );
|
|
|
100 |
$I = atan( sin( $rotI )
|
|
|
101 |
/ (cos( $this->b0 ) * cos( $rotI ) - sin( $this->b0 )
|
|
|
102 |
* tan( $rotB )) );
|
|
|
103 |
|
|
|
104 |
$lambda = $this->lambda0 + $I / $this->alpha;
|
|
|
105 |
|
|
|
106 |
$S = 0.0;
|
|
|
107 |
$phy = $b;
|
|
|
108 |
$prevPhy = -1000.0;
|
|
|
109 |
$iteration = 0;
|
|
|
110 |
while( abs( $phy - $prevPhy ) > 0.0000001 ) {
|
|
|
111 |
if( ++$iteration > 20 ) {
|
|
|
112 |
Proj4php::reportError( "omercFwdInfinity" );
|
|
|
113 |
return;
|
|
|
114 |
}
|
|
|
115 |
//S = log(tan(PI / 4.0 + phy / 2.0));
|
|
|
116 |
$S = 1.0
|
|
|
117 |
/ $this->alpha
|
|
|
118 |
* (log( tan( $PI / 4.0 + $b / 2.0 ) ) - $this->K)
|
|
|
119 |
+ $this->e
|
|
|
120 |
* log( tan( $PI / 4.0
|
|
|
121 |
+ asin( $this->e * sin( $phy ) )
|
|
|
122 |
/ 2.0 ) );
|
|
|
123 |
$prevPhy = $phy;
|
|
|
124 |
$phy = 2.0 * atan( exp( $S ) ) - $PI / 2.0;
|
|
|
125 |
}
|
|
|
126 |
|
|
|
127 |
$p->x = $lambda;
|
|
|
128 |
$p->y = $phy;
|
|
|
129 |
|
|
|
130 |
return $p;
|
|
|
131 |
}
|
|
|
132 |
|
|
|
133 |
}
|
|
|
134 |
|
|
|
135 |
Proj4php::$proj['somerc'] = new Proj4phpProjSomerc();
|