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